use crate::config::{AlignmentMode, ConsensusMode, PoaConfig};
use crate::error::PoaError;
use crate::types::{BubbleSite, Consensus, CoverageGap, GapKind, GraphStats};
use std::collections::HashMap;
#[cfg(test)]
use std::cell::Cell as StdCell;
#[cfg(test)]
thread_local! {
pub(crate) static SKIP_COUNTER: StdCell<u64> = StdCell::new(0);
pub(crate) static NODE_COUNTER: StdCell<u64> = StdCell::new(0);
}
#[cfg(test)]
pub(crate) fn reset_skip_counters() {
SKIP_COUNTER.with(|c| c.set(0));
NODE_COUNTER.with(|c| c.set(0));
}
#[cfg(test)]
pub(crate) fn skip_rate() -> f64 {
let skips = SKIP_COUNTER.with(|c| c.get());
let nodes = NODE_COUNTER.with(|c| c.get());
if nodes == 0 {
0.0
} else {
skips as f64 / nodes as f64
}
}
const UNSET: i32 = i32::MIN / 2;
const LOCK_EPS: usize = 5;
const LOCAL_POP_RADIUS: usize = 50;
const MINI_K: usize = 15;
const MINI_W: usize = 10;
const MINI_EPS_BASE: usize = 3;
const MINI_MIN_CHAIN: usize = 15;
struct AlignScratch {
lock_node_j: Vec<(usize, u32)>,
lock_exit_j: Vec<(usize, u32)>,
}
impl AlignScratch {
fn new() -> Self {
Self {
lock_node_j: Vec::new(),
lock_exit_j: Vec::new(),
}
}
fn clear(&mut self) {
self.lock_node_j.clear();
self.lock_exit_j.clear();
}
fn get_node_j(&self, node_idx: usize) -> Option<usize> {
self.lock_node_j
.binary_search_by_key(&node_idx, |&(idx, _)| idx)
.ok()
.map(|pos| self.lock_node_j[pos].1 as usize)
}
fn get_exit_j(&self, node_idx: usize) -> Option<usize> {
self.lock_exit_j
.binary_search_by_key(&node_idx, |&(idx, _)| idx)
.ok()
.map(|pos| self.lock_exit_j[pos].1 as usize)
}
}
#[inline]
fn edge_weight(nodes: &[Node], from: usize, to: usize) -> i32 {
nodes[from]
.out_edges
.iter()
.find(|&&(t, _)| t == to)
.map(|&(_, ew)| ew.total())
.unwrap_or(0)
}
#[inline]
fn safe_add(a: i32, b: i32) -> i32 {
if a == UNSET {
UNSET
} else {
a.saturating_add(b)
}
}
#[derive(Clone, Copy, Default, Debug, PartialEq, Eq)]
struct EdgeWeight {
matched: i32,
deleted: i32,
}
impl EdgeWeight {
fn total(&self) -> i32 {
self.matched + self.deleted
}
}
struct Node {
base: u8,
out_edges: Vec<(usize, EdgeWeight)>,
in_edges: Vec<usize>,
coverage: u32,
delete_count: u32,
nearest_fork: Option<(usize, usize)>,
}
pub struct PoaGraph {
nodes: Vec<Node>,
config: PoaConfig,
n_reads: usize,
reads: Vec<Vec<u8>>,
edge_reads: HashMap<(usize, usize), Vec<u32>>,
#[allow(dead_code)]
edge_delete_reads: HashMap<(usize, usize), Vec<u32>>,
bypass_edges: HashMap<usize, Vec<(usize, i32)>>,
warnings: usize,
cached_spine: Vec<(usize, u8, i32)>,
spine_updated_at: usize,
spine_interval: usize,
align_scratch: AlignScratch,
spine_mers: HashMap<u64, u32>,
fork_arm_index: HashMap<usize, HashMap<Vec<u8>, usize>>,
used_band_retry: bool,
}
const VIRTUAL: u32 = u32::MAX;
#[derive(Clone, Copy, PartialEq, Eq)]
enum State {
M,
I,
D,
}
#[derive(Clone, Copy)]
struct Cell {
score: i32,
pred_t: u32,
}
impl Cell {
fn unset() -> Self {
Cell {
score: UNSET,
pred_t: 0,
}
}
}
#[derive(Clone, Copy, Debug)]
pub enum AlignOp {
Match(usize),
Insert(u8),
Delete(usize),
}
fn push_node(nodes: &mut Vec<Node>, base: u8) -> usize {
let idx = nodes.len();
nodes.push(Node {
base,
out_edges: Vec::new(),
in_edges: Vec::new(),
coverage: 0,
delete_count: 0,
nearest_fork: None,
});
idx
}
const FORK_PROPAGATE_MAX_DEPTH: usize = ARM_MAX_DEPTH;
fn propagate_fork_if_new(nodes: &mut [Node], from: usize) {
if nodes[from].out_edges.len() != 2 {
return;
}
let orig_target = nodes[from].out_edges[0].0;
let mut cur = orig_target;
for _ in 0..FORK_PROPAGATE_MAX_DEPTH {
if nodes[cur].in_edges.len() > 1 {
break;
}
nodes[cur].nearest_fork = Some((from, orig_target));
match nodes[cur].out_edges.as_slice() {
[(next, _)] => cur = *next,
_ => break,
}
}
}
fn set_new_node_own_fork(nodes: &mut [Node], p: usize, new_idx: usize) {
if nodes[p].out_edges.len() >= 2 {
nodes[new_idx].nearest_fork = Some((p, new_idx));
} else {
nodes[new_idx].nearest_fork = nodes[p].nearest_fork;
}
}
fn add_edge(nodes: &mut [Node], from: usize, to: usize) {
nodes[from].out_edges.push((
to,
EdgeWeight {
matched: 1,
deleted: 0,
},
));
nodes[to].in_edges.push(from);
propagate_fork_if_new(nodes, from);
}
#[must_use]
fn increment_or_add_edge(nodes: &mut [Node], from: usize, to: usize, is_delete: bool) -> bool {
debug_assert_ne!(
from, to,
"increment_or_add_edge: refusing to create/increment a literal self-loop"
);
for (succ, ew) in nodes[from].out_edges.iter_mut() {
if *succ == to {
if is_delete {
ew.deleted += 1;
} else {
ew.matched += 1;
}
return false;
}
}
let ew = if is_delete {
EdgeWeight {
matched: 0,
deleted: 1,
}
} else {
EdgeWeight {
matched: 1,
deleted: 0,
}
};
nodes[from].out_edges.push((to, ew));
nodes[to].in_edges.push(from);
true
}
fn topological_order(nodes: &[Node]) -> (Vec<usize>, Vec<usize>) {
let n = nodes.len();
let mut in_deg: Vec<usize> = nodes.iter().map(|nd| nd.in_edges.len()).collect();
let mut queue: std::collections::VecDeque<usize> = (0..n).filter(|&i| in_deg[i] == 0).collect();
let mut topo: Vec<usize> = Vec::with_capacity(n);
while let Some(u) = queue.pop_front() {
topo.push(u);
for &(v, _) in &nodes[u].out_edges {
in_deg[v] -= 1;
if in_deg[v] == 0 {
queue.push_back(v);
}
}
}
let mut rank_of = vec![0usize; n];
for (t, &node_idx) in topo.iter().enumerate() {
rank_of[node_idx] = t;
}
(topo, rank_of)
}
fn materialize_arm_len(nodes: &[Node], start: usize, max_depth: usize) -> usize {
let mut len = 0;
let mut cur = start;
for _ in 0..max_depth {
if len >= 1 && nodes[cur].in_edges.len() > 1 {
break;
}
len += 1;
match nodes[cur].out_edges.as_slice() {
[(next, _)] => cur = *next,
_ => break,
}
}
len
}
fn real_in_edge_count(nodes: &[Node], node: usize, sig_threshold: i32) -> usize {
nodes[node]
.in_edges
.iter()
.filter(|&&p| {
nodes[p]
.out_edges
.iter()
.find(|&&(to, _)| to == node)
.is_some_and(|&(_, ew)| ew.matched >= sig_threshold)
})
.count()
}
fn materialize_arm_len_tolerant(
nodes: &[Node],
start: usize,
max_depth: usize,
sig_threshold: i32,
) -> usize {
let mut len = 0;
let mut cur = start;
for _ in 0..max_depth {
if len >= 1 && real_in_edge_count(nodes, cur, sig_threshold) > 1 {
break;
}
len += 1;
match nodes[cur].out_edges.as_slice() {
[] => break,
[(next, _)] => cur = *next,
edges => {
let real_arms = edges
.iter()
.filter(|&&(_, ew)| ew.matched >= sig_threshold)
.count();
if real_arms > 1 {
break; }
let (next, _) = edges.iter().max_by_key(|&&(_, ew)| ew.matched).unwrap();
cur = *next;
}
}
}
len
}
fn compute_bubble_ranges(nodes: &[Node], topo: &[usize]) -> Vec<Option<(usize, usize)>> {
let n = topo.len();
let nn = nodes.len();
let mut ranges: Vec<Option<(usize, usize)>> = vec![None; n];
let mut gen_mark = vec![0u32; nn];
let mut current_gen = 0u32;
let mut t = 0;
while t < n {
if nodes[topo[t]].out_edges.len() >= 2 {
current_gen = current_gen.wrapping_add(1);
let entry_node = topo[t];
gen_mark[entry_node] = current_gen;
let mut outstanding = nodes[entry_node].out_edges.len();
let mut exit_t = None;
for (tt, &node_idx) in topo.iter().enumerate().skip(t + 1) {
let bubble_in = nodes[node_idx]
.in_edges
.iter()
.filter(|&&p| gen_mark[p] == current_gen)
.count();
if bubble_in == 0 {
continue;
}
gen_mark[node_idx] = current_gen;
outstanding -= bubble_in;
if outstanding == 0 {
exit_t = Some(tt);
break;
}
outstanding += nodes[node_idx].out_edges.len();
}
if let Some(et) = exit_t {
ranges[t..=et].fill(Some((t, et)));
t = et + 1;
} else {
ranges[t..n].fill(Some((t, n - 1)));
break;
}
} else {
t += 1;
}
}
ranges
}
fn spine_diff(old: &[(usize, u8, i32)], new: &[(usize, u8, i32)]) -> usize {
if old.is_empty() {
return usize::MAX;
}
let len_diff = old.len().abs_diff(new.len());
let base_diffs = old[..old.len().min(new.len())]
.iter()
.zip(new.iter())
.filter(|(o, n)| o.1 != n.1)
.count();
base_diffs + len_diff
}
const SPINE_STABLE_THRESHOLD: usize = 3;
const SPINE_MAX_INTERVAL: usize = 32;
const LOOKAHEAD_K: usize = 5;
const LOOKAHEAD_MARGIN: i32 = 2;
const ARM_MAX_DEPTH: usize = 4096;
fn collect_arm_nodes(
nodes: &[Node],
rank_of: &[usize],
start_idx: usize,
exit_t: usize,
out: &mut Vec<usize>,
) -> bool {
out.clear();
let mut cur = start_idx;
for _ in 0..ARM_MAX_DEPTH {
if rank_of[cur] == exit_t {
return true;
}
out.push(cur);
match nodes[cur].out_edges.as_slice() {
[(next, _)] => cur = *next,
_ => return false,
}
}
false
}
fn score_arm_prefix(
arm: &[usize],
nodes: &[Node],
query: &[u8],
j: usize,
match_score: i32,
mismatch_score: i32,
) -> i32 {
let len = arm
.len()
.min(LOOKAHEAD_K)
.min(query.len().saturating_sub(j));
arm[..len]
.iter()
.enumerate()
.map(|(i, &nidx)| {
if nodes[nidx].base == query[j + i] {
match_score
} else {
mismatch_score
}
})
.sum()
}
const SLIDE_MIN_CONSEC: usize = 2;
fn slide_lock(
all_arms: &[Vec<usize>],
nodes: &[Node],
query: &[u8],
j_entry: usize,
) -> Option<usize> {
let n_arms = all_arms.len();
let remaining = query.len().saturating_sub(j_entry);
if remaining == 0 || n_arms < 2 {
return None;
}
let min_arm_len = all_arms.iter().map(|a| a.len()).min().unwrap_or(0);
if min_arm_len < LOOKAHEAD_K {
return None;
}
let max_steps = all_arms
.iter()
.map(|a| a.len())
.max()
.unwrap_or(0)
.min(remaining);
let mut alive = vec![true; n_arms];
let mut alive_count = n_arms;
let mut consec_unique = 0usize;
let mut consec_winner: Option<usize> = None;
for i in 0..max_steps {
let any_continuing = (0..n_arms).any(|idx| alive[idx] && all_arms[idx].len() > i);
if any_continuing {
for idx in 0..n_arms {
if alive[idx] && all_arms[idx].len() <= i {
alive[idx] = false;
alive_count -= 1;
}
}
}
if alive_count <= 1 {
break;
}
let q = query[j_entry + i];
let matched: Vec<bool> = (0..n_arms)
.map(|idx| alive[idx] && all_arms[idx].len() > i && nodes[all_arms[idx][i]].base == q)
.collect();
let match_count = matched.iter().filter(|&&m| m).count();
let mismatch_count = (0..n_arms)
.filter(|&idx| alive[idx] && !matched[idx])
.count();
if match_count == 1 && mismatch_count > 0 {
let candidate = matched.iter().position(|&m| m).unwrap();
if consec_winner == Some(candidate) {
consec_unique += 1;
} else {
consec_unique = 1;
consec_winner = Some(candidate);
}
if consec_unique >= SLIDE_MIN_CONSEC {
for idx in 0..n_arms {
if alive[idx] && !matched[idx] {
alive[idx] = false;
alive_count -= 1;
}
}
break;
}
} else {
consec_unique = 0;
consec_winner = None;
}
}
if alive_count == 1 {
alive.iter().position(|&a| a)
} else {
None
}
}
#[inline]
fn encode_base(b: u8) -> Option<u64> {
match b {
b'A' | b'a' => Some(0),
b'C' | b'c' => Some(1),
b'G' | b'g' => Some(2),
b'T' | b't' => Some(3),
_ => None,
}
}
fn compute_minimizers(seq: &[u8], k: usize, w: usize) -> Vec<(u64, usize)> {
let n = seq.len();
if n < k {
return vec![];
}
let n_kmers = n - k + 1;
let mask = if k * 2 < 64 {
(1u64 << (k * 2)) - 1
} else {
u64::MAX
};
let mut kmer_hashes: Vec<Option<u64>> = Vec::with_capacity(n_kmers);
let mut hash: u64 = 0;
let mut valid: usize = 0;
for (i, &b) in seq.iter().enumerate() {
match encode_base(b) {
Some(bits) => {
hash = ((hash << 2) | bits) & mask;
valid += 1;
}
None => {
hash = 0;
valid = 0;
}
}
if i + 1 >= k {
kmer_hashes.push(if valid >= k { Some(hash) } else { None });
}
}
let win = w.min(n_kmers);
let mut result: Vec<(u64, usize)> = Vec::new();
let mut last: Option<(u64, usize)> = None;
for start in 0..=(n_kmers - win) {
let mut best_hash = u64::MAX;
let mut best_pos = start;
for (off, entry) in kmer_hashes[start..start + win].iter().enumerate() {
if let Some(h) = entry {
if *h < best_hash {
best_hash = *h;
best_pos = start + off;
}
}
}
if best_hash == u64::MAX {
continue;
}
let entry = (best_hash, best_pos);
if last != Some(entry) {
last = Some(entry);
result.push(entry);
}
}
result
}
fn build_spine_mers(spine: &[(usize, u8, i32)], k: usize, w: usize) -> HashMap<u64, u32> {
let spine_seq: Vec<u8> = spine.iter().map(|&(_, b, _)| b).collect();
let mers = compute_minimizers(&spine_seq, k, w);
let mut counts: HashMap<u64, (u32, u32)> = HashMap::new();
for (hash, pos) in mers {
let e = counts.entry(hash).or_insert((0, pos as u32));
e.0 += 1;
}
counts
.into_iter()
.filter(|(_, (c, _))| *c == 1)
.map(|(h, (_, p))| (h, p))
.collect()
}
fn lis_indices(vals: &[usize]) -> Vec<usize> {
let n = vals.len();
if n == 0 {
return vec![];
}
let mut tails: Vec<usize> = Vec::new();
let mut tail_idx: Vec<usize> = Vec::new();
let mut pred: Vec<Option<usize>> = vec![None; n];
for i in 0..n {
let v = vals[i];
let p = tails.partition_point(|&t| t < v);
if p == tails.len() {
tails.push(v);
tail_idx.push(i);
} else {
tails[p] = v;
tail_idx[p] = i;
}
pred[i] = if p > 0 { Some(tail_idx[p - 1]) } else { None };
}
let lis_len = tails.len();
let mut chain = Vec::with_capacity(lis_len);
let mut cur = Some(tail_idx[lis_len - 1]);
while let Some(i) = cur {
chain.push(i);
cur = pred[i];
}
chain.reverse();
chain
}
fn build_anchor_chain(
read_mers: &[(u64, usize)],
spine_mers: &HashMap<u64, u32>,
spine: &[(usize, u8, i32)],
rank_of: &[usize],
) -> Vec<(usize, usize)> {
let mut candidates: Vec<(usize, usize)> = read_mers
.iter()
.filter_map(|&(hash, read_pos)| {
spine_mers.get(&hash).map(|&sr| {
let node_idx = spine[sr as usize].0;
(read_pos, rank_of[node_idx])
})
})
.collect();
if candidates.is_empty() {
return vec![];
}
candidates.sort_unstable_by_key(|&(_, t)| t);
let read_positions: Vec<usize> = candidates.iter().map(|&(r, _)| r).collect();
let idxs = lis_indices(&read_positions);
idxs.iter().map(|&i| candidates[i]).collect()
}
fn anchor_j_bounds(
t: usize,
anchors: &[(usize, usize)], l: usize,
) -> Option<(usize, usize)> {
if anchors.is_empty() {
return None;
}
let pos = anchors.partition_point(|&(_, at)| at <= t);
if pos > 0 && anchors[pos - 1].1 == t {
let r = anchors[pos - 1].0;
let lo = r.saturating_sub(MINI_EPS_BASE).max(1);
let hi = (r + MINI_EPS_BASE).min(l);
Some((lo, hi))
} else {
None
}
}
#[allow(clippy::too_many_arguments)]
#[inline]
fn anchor_refine_spine(
j_lo: usize,
j_hi: usize,
t: usize,
anchors: &[(usize, usize)],
l: usize,
is_source: bool,
is_spine: bool,
in_bubble: bool,
is_locked: bool,
j_center: usize,
) -> (usize, usize) {
if anchors.is_empty() || is_source || !is_spine || in_bubble || is_locked {
return (j_lo, j_hi);
}
match anchor_j_bounds(t, anchors, l) {
Some((alo, ahi)) => {
if ahi < j_center {
return (j_lo, j_hi);
}
let nlo = j_lo.max(alo);
let nhi = j_hi.min(ahi);
if nlo <= nhi { (nlo, nhi) } else { (j_lo, j_hi) }
}
None => (j_lo, j_hi),
}
}
#[inline(always)]
fn gs(mat: &[Cell], t: usize, j: usize, j_lo: usize, j_hi: usize, rw: usize) -> i32 {
if j == 0 || j < j_lo || j > j_hi {
UNSET
} else {
mat[t * rw + (j - j_lo)].score
}
}
#[inline(always)]
fn gsd(
del: &[Cell],
del0: &[Cell],
t: usize,
j: usize,
j_lo: usize,
j_hi: usize,
rw: usize,
) -> i32 {
if j == 0 {
del0[t].score
} else if j < j_lo || j > j_hi {
UNSET
} else {
del[t * rw + (j - j_lo)].score
}
}
#[allow(clippy::too_many_arguments)]
fn align(
nodes: &[Node],
topo: &[usize],
rank_of: &[usize],
spine: &[(usize, u8, i32)],
query: &[u8],
cfg: &PoaConfig,
scratch: &mut AlignScratch,
anchors: &[(usize, usize)], ) -> Result<Vec<AlignOp>, PoaError> {
scratch.clear();
let n = topo.len();
let l = query.len();
if cfg.warn_on_long_unbanded && cfg.band_width == 0 && !cfg.adaptive_band && l > 1000 {
eprintln!(
"poa-consensus: warning: unbanded alignment (band_width=0, \
adaptive_band=false) on a {l} bp read -- this scales as \
O(read_len * graph_len); consider a banded or adaptive PoaConfig \
for large graphs, or set warn_on_long_unbanded=false to suppress"
);
}
let go = cfg.gap_open;
let ge = cfg.gap_extend;
let semi = cfg.alignment_mode == AlignmentMode::SemiGlobal;
let nn = nodes.len();
let mut on_spine = vec![false; nn];
let mut spine_prev: Vec<Option<usize>> = vec![None; nn];
for (sp, &(node_idx, _, _)) in spine.iter().enumerate() {
on_spine[node_idx] = true;
if sp > 0 {
spine_prev[node_idx] = Some(spine[sp - 1].0);
}
}
let bubble_ranges = compute_bubble_ranges(nodes, topo);
let diag_skip_disabled = !cfg.multi_allele;
let mut bubble_entry_j = vec![0usize; n];
let mut best_j_per_t = vec![0usize; n];
let mut lookahead_skip = vec![false; nn];
let spine_margin: usize = if cfg.adaptive_band {
let w = cfg.adaptive_band_b + (cfg.adaptive_band_f * l as f32).ceil() as usize;
let w = if cfg.band_width > 0 {
w.max(cfg.band_width)
} else {
w
};
w.max(4)
} else if cfg.band_width > 0 {
cfg.band_width
} else {
l
};
let row_width = (2 * spine_margin + 2).min(l).max(1);
let mut del0 = vec![Cell::unset(); n];
let bsize = n * row_width;
let mut m = vec![Cell::unset(); bsize];
let mut ins = vec![Cell::unset(); bsize];
let mut del = vec![Cell::unset(); bsize];
let mut j_lo_arr = vec![1usize; n];
let mut j_hi_arr = vec![0usize; n];
let mut arm_scratch: Vec<usize> = Vec::new();
let mut all_arms: Vec<Vec<usize>> = Vec::new();
for (t, &node_idx) in topo.iter().enumerate() {
if lookahead_skip[node_idx] {
continue;
}
let node_base = nodes[node_idx].base;
let is_source = nodes[node_idx].in_edges.is_empty();
let j_pred_max = if is_source {
0
} else {
nodes[node_idx]
.in_edges
.iter()
.map(|&p| best_j_per_t[rank_of[p]])
.max()
.unwrap_or(0)
};
let (j_lo, j_hi) = if is_source {
(1usize, row_width.min(l))
} else if let Some(j_center) = scratch.get_node_j(node_idx) {
let lo = j_center.saturating_sub(LOCK_EPS).max(1);
let hi = (j_center + LOCK_EPS).min(l);
(lo, hi)
} else if let Some(j_center) = scratch.get_exit_j(node_idx) {
let lo = j_center.saturating_sub(spine_margin).max(1);
let hi = (j_center + spine_margin).min(l);
(lo, hi)
} else {
match bubble_ranges[t] {
Some((entry_t, exit_t)) => {
if t == entry_t {
bubble_entry_j[entry_t] = j_pred_max;
}
let bej = bubble_entry_j[entry_t];
let bubble_span = exit_t.saturating_sub(entry_t) + 1;
let lo = bej.saturating_sub(spine_margin).max(1);
let hi = (bej + bubble_span.min(spine_margin) + spine_margin).min(l);
(lo, hi)
}
None => {
let j_center = j_pred_max.saturating_add(1);
let lo = j_center.saturating_sub(spine_margin).max(1);
let hi = (j_center + spine_margin).min(l);
(lo, hi)
}
}
};
let j_hi = j_hi.min(j_lo + row_width - 1);
let (j_lo, j_hi) = anchor_refine_spine(
j_lo,
j_hi,
t,
anchors,
l,
is_source,
on_spine[node_idx],
bubble_ranges[t].is_some(),
scratch.get_node_j(node_idx).is_some() || scratch.get_exit_j(node_idx).is_some(),
j_pred_max.saturating_add(1),
);
j_lo_arr[t] = j_lo;
j_hi_arr[t] = j_hi;
#[cfg(test)]
if !is_source {
NODE_COUNTER.with(|c| c.set(c.get() + 1));
}
if !is_source && on_spine[node_idx] && !diag_skip_disabled {
if let Some(prev_sp) = spine_prev[node_idx] {
let active_pred_ok = match nodes[node_idx].in_edges.as_slice() {
[p] => *p == prev_sp,
_ => {
nodes[node_idx]
.in_edges
.iter()
.all(|&p| p == prev_sp || lookahead_skip[p])
&& nodes[node_idx].in_edges.contains(&prev_sp)
}
};
if active_pred_ok {
let pt = rank_of[prev_sp];
let bj = best_j_per_t[pt];
if bj < l && bj + 1 >= j_lo && bj < j_hi && node_base == query[bj] {
let m_prev = gs(&m, pt, bj, j_lo_arr[pt], j_hi_arr[pt], row_width);
let i_prev = gs(&ins, pt, bj, j_lo_arr[pt], j_hi_arr[pt], row_width);
let d_prev =
gsd(&del, &del0, pt, bj, j_lo_arr[pt], j_hi_arr[pt], row_width);
let pred_is_source = nodes[topo[pt]].in_edges.is_empty();
let i_ok = i_prev == UNSET || (pred_is_source && m_prev > i_prev);
let d_ok = d_prev == UNSET || (pred_is_source && m_prev > d_prev);
if m_prev != UNSET && i_ok && d_ok {
let do_skip = match nodes[node_idx].out_edges.as_slice() {
[_] => true,
_ if bj + 1 < l => {
let next_q = query[bj + 1];
let mut spine_succ = None;
let mut resolved = true;
for &(s, _) in &nodes[node_idx].out_edges {
if on_spine[s] && spine_prev[s] == Some(node_idx) {
if nodes[s].base == next_q {
spine_succ = Some(s);
} else {
resolved = false;
break;
}
} else if !lookahead_skip[s] && nodes[s].base == next_q {
resolved = false;
break;
}
}
if resolved {
if let Some(ss) = spine_succ {
for &(s, _) in &nodes[node_idx].out_edges {
if s != ss && !lookahead_skip[s] {
lookahead_skip[s] = true;
let mut cur = s;
for _ in 0..ARM_MAX_DEPTH {
match nodes[cur].out_edges.as_slice() {
[(next, _)]
if !on_spine[*next]
&& !lookahead_skip[*next] =>
{
lookahead_skip[*next] = true;
cur = *next;
}
_ => break,
}
}
}
}
true
} else {
false
}
} else {
false
}
}
_ => false,
};
if do_skip {
let score = m_prev + cfg.match_score;
m[t * row_width + (bj + 1 - j_lo)] = Cell {
score,
pred_t: pt as u32,
};
best_j_per_t[t] = bj + 1;
#[cfg(test)]
SKIP_COUNTER.with(|c| c.set(c.get() + 1));
continue;
}
}
}
}
}
}
{
let (mut best, mut best_pred) = (UNSET, 0u32);
if is_source {
let val = go + ge;
if val > best {
best = val;
best_pred = VIRTUAL;
}
}
for &p in &nodes[node_idx].in_edges {
let pt = rank_of[p];
let vd = safe_add(del0[pt].score, ge);
if vd > best {
best = vd;
best_pred = pt as u32;
}
}
if best != UNSET {
del0[t] = Cell {
score: best,
pred_t: best_pred,
};
}
}
let mut row_best_j = 0usize;
let mut row_best_score = UNSET;
for j in j_lo..=j_hi {
let q_base = query[j - 1];
let sc = if node_base == q_base {
cfg.match_score
} else {
cfg.mismatch_score
};
let ixcur = t * row_width + (j - j_lo);
{
let (mut best, mut best_pred) = (UNSET, VIRTUAL);
if j == 1 && (is_source || semi) && sc > best {
best = sc;
best_pred = VIRTUAL;
}
if is_source && j > 1 {
let val = safe_add(go + (j as i32 - 1) * ge, sc);
if val != UNSET && val > best {
best = val;
best_pred = VIRTUAL;
}
}
let mut best_edge_w = i32::MAX;
for &p in &nodes[node_idx].in_edges {
let pt = rank_of[p];
let ew = edge_weight(nodes, p, node_idx);
let vm = safe_add(gs(&m, pt, j - 1, j_lo_arr[pt], j_hi_arr[pt], row_width), sc);
if vm != UNSET && (vm > best || (vm == best && ew > best_edge_w)) {
best = vm;
best_pred = pt as u32;
best_edge_w = ew;
}
let vi = safe_add(
gs(&ins, pt, j - 1, j_lo_arr[pt], j_hi_arr[pt], row_width),
sc,
);
if vi != UNSET && (vi > best || (vi == best && ew > best_edge_w)) {
best = vi;
best_pred = pt as u32;
best_edge_w = ew;
}
let vd = safe_add(
gsd(
&del,
&del0,
pt,
j - 1,
j_lo_arr[pt],
j_hi_arr[pt],
row_width,
),
sc,
);
if vd != UNSET && (vd > best || (vd == best && ew > best_edge_w)) {
best = vd;
best_pred = pt as u32;
best_edge_w = ew;
}
}
if best != UNSET {
m[ixcur] = Cell {
score: best,
pred_t: best_pred,
};
if best > row_best_score {
row_best_score = best;
row_best_j = j;
}
}
}
{
let (mut best, mut best_pred) = (UNSET, VIRTUAL);
if is_source && j == 1 {
let val = go + ge;
if val > best {
best = val;
best_pred = VIRTUAL;
}
}
if j > j_lo {
let ixprev = t * row_width + (j - 1 - j_lo);
let vm = safe_add(m[ixprev].score, go + ge);
if vm != UNSET && vm > best {
best = vm;
best_pred = t as u32;
}
let vi = safe_add(ins[ixprev].score, ge);
if vi != UNSET && vi > best {
best = vi;
best_pred = t as u32;
}
}
if best != UNSET {
ins[ixcur] = Cell {
score: best,
pred_t: best_pred,
};
}
}
{
let (mut best, mut best_pred) = (UNSET, 0u32);
let mut best_edge_w = i32::MAX;
for &p in &nodes[node_idx].in_edges {
let pt = rank_of[p];
let ew = edge_weight(nodes, p, node_idx);
let vm = safe_add(
gs(&m, pt, j, j_lo_arr[pt], j_hi_arr[pt], row_width),
go + ge,
);
if vm != UNSET && (vm > best || (vm == best && ew > best_edge_w)) {
best = vm;
best_pred = pt as u32;
best_edge_w = ew;
}
let vi = safe_add(
gs(&ins, pt, j, j_lo_arr[pt], j_hi_arr[pt], row_width),
go + ge,
);
if vi != UNSET && (vi > best || (vi == best && ew > best_edge_w)) {
best = vi;
best_pred = pt as u32;
best_edge_w = ew;
}
let vd = safe_add(
gsd(&del, &del0, pt, j, j_lo_arr[pt], j_hi_arr[pt], row_width),
ge,
);
if vd != UNSET && (vd > best || (vd == best && ew > best_edge_w)) {
best = vd;
best_pred = pt as u32;
best_edge_w = ew;
}
}
if best != UNSET {
del[ixcur] = Cell {
score: best,
pred_t: best_pred,
};
}
}
}
if row_best_score != UNSET {
best_j_per_t[t] = row_best_j;
} else {
best_j_per_t[t] = if is_source {
1
} else {
j_pred_max.saturating_add(1).min(l)
};
}
if let Some((entry_t, exit_t)) = bubble_ranges[t] {
if t == entry_t && nodes[node_idx].out_edges.len() >= 2 {
let j_entry = best_j_per_t[t];
if j_entry < l {
all_arms.clear();
let mut complex = false;
for &(arm_start, _) in &nodes[node_idx].out_edges {
if collect_arm_nodes(nodes, rank_of, arm_start, exit_t, &mut arm_scratch) {
all_arms.push(arm_scratch.clone());
} else {
complex = true;
break;
}
}
if !complex && all_arms.len() >= 2 {
let min_scored = all_arms
.iter()
.map(|arm| {
arm.len()
.min(LOOKAHEAD_K)
.min(query.len().saturating_sub(j_entry))
})
.min()
.unwrap_or(0);
let winner = if min_scored >= LOOKAHEAD_K {
let scores: Vec<i32> = all_arms
.iter()
.map(|arm| {
score_arm_prefix(
arm,
nodes,
query,
j_entry,
cfg.match_score,
cfg.mismatch_score,
)
})
.collect();
let best = *scores.iter().max().unwrap();
let winner_count = scores.iter().filter(|&&s| s == best).count();
if winner_count == 1 {
let second_best = scores
.iter()
.copied()
.filter(|&s| s != best)
.max()
.unwrap_or(i32::MIN);
if best - second_best >= LOOKAHEAD_MARGIN {
scores.iter().position(|&s| s == best)
} else {
None
}
} else {
None
}
} else {
None
};
let winner =
winner.or_else(|| slide_lock(&all_arms, nodes, query, j_entry));
if let Some(w) = winner {
for (arm_idx, arm) in all_arms.iter().enumerate() {
if arm_idx != w {
for &losing_node in arm {
lookahead_skip[losing_node] = true;
}
}
}
let arm = &all_arms[w];
let arm_len = arm.len();
for (d, &win_node) in arm.iter().enumerate() {
let j_c = (j_entry + d + 1).min(l) as u32;
scratch.lock_node_j.push((win_node, j_c));
}
scratch.lock_node_j.sort_unstable_by_key(|&(idx, _)| idx);
let exit_node = topo[exit_t];
let exit_j = (j_entry + arm_len).min(l) as u32;
match scratch
.lock_exit_j
.binary_search_by_key(&exit_node, |&(idx, _)| idx)
{
Ok(pos) => {
scratch.lock_exit_j[pos].1 =
exit_j.max(scratch.lock_exit_j[pos].1)
}
Err(pos) => scratch.lock_exit_j.insert(pos, (exit_node, exit_j)),
}
}
}
}
}
}
}
let terminal_best = (0..n)
.flat_map(|t| {
let jlo = j_lo_arr[t];
let jhi = j_hi_arr[t];
let sm = gs(&m, t, l, jlo, jhi, row_width);
let si = gs(&ins, t, l, jlo, jhi, row_width);
let sd = gsd(&del, &del0, t, l, jlo, jhi, row_width);
let best_sc = [sm, si, sd].into_iter().filter(|&s| s != UNSET).max();
best_sc.map(|sc| {
let st = if sm == sc {
State::M
} else if si == sc {
State::I
} else {
State::D
};
(t, st, sc)
})
})
.max_by_key(|&(_, _, sc)| sc)
.map(|(t, s, _)| (t, s));
let (best_t, best_state) = match terminal_best {
Some(result) => result,
None => {
let max_best_j = best_j_per_t.iter().copied().max().unwrap_or(0);
let shortfall = l.saturating_sub(max_best_j);
let required = spine_margin.saturating_add(shortfall);
return Err(PoaError::BandTooNarrow {
configured: spine_margin,
required,
});
}
};
let mut ops: Vec<AlignOp> = Vec::with_capacity(l + n / 4);
let mut t = best_t;
let mut j = l;
let mut cur_state = best_state;
let max_ops = l.saturating_add(n).saturating_add(16);
loop {
if ops.len() > max_ops {
return Err(PoaError::BandTooNarrow {
configured: spine_margin,
required: spine_margin.saturating_mul(2).max(l),
});
}
let cell = {
let jlo = j_lo_arr[t];
let jhi = j_hi_arr[t];
match cur_state {
State::M => {
if j < jlo || j > jhi {
Cell::unset()
} else {
m[t * row_width + (j - jlo)]
}
}
State::I => {
if j < jlo || j > jhi {
Cell::unset()
} else {
ins[t * row_width + (j - jlo)]
}
}
State::D => {
if j == 0 {
del0[t]
} else if j < jlo || j > jhi {
Cell::unset()
} else {
del[t * row_width + (j - jlo)]
}
}
}
};
if cell.score == UNSET {
break;
}
match cur_state {
State::M => {
ops.push(AlignOp::Match(topo[t]));
if cell.pred_t == VIRTUAL {
for k in (1..j).rev() {
ops.push(AlignOp::Insert(query[k - 1]));
}
break;
}
t = cell.pred_t as usize;
j -= 1;
cur_state = best_prev_state_banded(
&m, &ins, &del, &del0, t, j, &j_lo_arr, &j_hi_arr, row_width,
);
}
State::I => {
ops.push(AlignOp::Insert(query[j - 1]));
j -= 1;
if cell.pred_t == VIRTUAL {
for k in (1..j).rev() {
ops.push(AlignOp::Insert(query[k - 1]));
}
break;
}
let sm = gs(&m, t, j, j_lo_arr[t], j_hi_arr[t], row_width);
let si = gs(&ins, t, j, j_lo_arr[t], j_hi_arr[t], row_width);
cur_state = if sm >= si { State::M } else { State::I };
}
State::D => {
ops.push(AlignOp::Delete(topo[t]));
if cell.pred_t == VIRTUAL {
break;
}
t = cell.pred_t as usize;
cur_state = best_prev_state_banded(
&m, &ins, &del, &del0, t, j, &j_lo_arr, &j_hi_arr, row_width,
);
}
}
if j == 0 && cur_state != State::D {
break;
}
}
ops.reverse();
Ok(ops)
}
#[allow(clippy::too_many_arguments)]
fn align_with_retry(
nodes: &[Node],
topo: &[usize],
rank_of: &[usize],
spine: &[(usize, u8, i32)],
query: &[u8],
cfg: &PoaConfig,
scratch: &mut AlignScratch,
anchors: &[(usize, usize)],
) -> Result<(Vec<AlignOp>, bool), PoaError> {
match align(nodes, topo, rank_of, spine, query, cfg, scratch, anchors) {
Ok(ops) => Ok((ops, false)),
Err(PoaError::BandTooNarrow { required, .. }) => {
let mut cfg2 = cfg.clone();
cfg2.band_width = required;
cfg2.adaptive_band = false;
match align(nodes, topo, rank_of, spine, query, &cfg2, scratch, &[]) {
Ok(ops) => Ok((ops, true)),
Err(PoaError::BandTooNarrow { .. }) => {
let mut cfg3 = cfg.clone();
cfg3.band_width = 0;
cfg3.adaptive_band = false;
align(nodes, topo, rank_of, spine, query, &cfg3, scratch, &[])
.map(|ops| (ops, true))
}
Err(other) => Err(other),
}
}
Err(other) => Err(other),
}
}
#[allow(clippy::too_many_arguments)]
#[inline]
fn best_prev_state_banded(
m: &[Cell],
ins: &[Cell],
del: &[Cell],
del0: &[Cell],
t: usize,
j: usize,
j_lo_arr: &[usize],
j_hi_arr: &[usize],
row_width: usize,
) -> State {
let jlo = j_lo_arr[t];
let jhi = j_hi_arr[t];
let sm = gs(m, t, j, jlo, jhi, row_width);
let si = gs(ins, t, j, jlo, jhi, row_width);
let sd = gsd(del, del0, t, j, jlo, jhi, row_width);
if sm != UNSET && sm >= si && sm >= sd {
State::M
} else if si != UNSET && si >= sd {
State::I
} else {
State::D
}
}
fn verify_reuse_chain(nodes: &[Node], edit: &[u8], start: usize) -> Option<Vec<usize>> {
debug_assert!(!edit.is_empty());
let mut chain = Vec::with_capacity(edit.len());
let mut cur = start;
if nodes[cur].base != edit[0] {
return None;
}
chain.push(cur);
for &b in &edit[1..] {
match nodes[cur].out_edges.as_slice() {
[(next, _)] if nodes[*next].base == b => {
cur = *next;
chain.push(cur);
}
_ => return None,
}
}
Some(chain)
}
fn reuse_would_collide(chain: &[usize], rest: &[AlignOp]) -> bool {
rest.iter().any(|op| match *op {
AlignOp::Match(idx) | AlignOp::Delete(idx) => chain.contains(&idx),
AlignOp::Insert(_) => false,
})
}
fn reuse_would_create_back_edge(rank_of: &[usize], chain: &[usize], rest: &[AlignOp]) -> bool {
let Some(&last_rank) = chain.last().and_then(|&idx| rank_of.get(idx)) else {
return false;
};
for op in rest {
if let AlignOp::Match(idx) | AlignOp::Delete(idx) = *op {
return match rank_of.get(idx) {
Some(&r) => last_rank >= r,
None => false,
};
}
}
false
}
fn commit_reuse_chain(
nodes: &mut [Node],
edge_reads: &mut HashMap<(usize, usize), Vec<u32>>,
fork: usize,
chain: &[usize],
read_idx: u32,
) -> usize {
let mut prev_node = fork;
for &node in chain {
let created = increment_or_add_edge(nodes, prev_node, node, false);
debug_assert!(
!created,
"reuse chain traversed an edge that did not already exist"
);
edge_reads
.entry((prev_node, node))
.or_default()
.push(read_idx);
nodes[node].coverage += 1;
prev_node = node;
}
*chain.last().unwrap()
}
#[allow(clippy::too_many_arguments)]
fn try_reuse_arm(
nodes: &mut [Node],
edge_reads: &mut HashMap<(usize, usize), Vec<u32>>,
fork_arm_index: &HashMap<usize, HashMap<Vec<u8>, usize>>,
rank_of: &[usize],
fork: usize,
edit: &[u8],
rest_ops: &[AlignOp],
read_idx: u32,
) -> Option<usize> {
let start = *fork_arm_index.get(&fork)?.get(edit)?;
if fork >= rank_of.len() || start >= rank_of.len() || rank_of[fork] >= rank_of[start] {
return None;
}
let chain = verify_reuse_chain(nodes, edit, start)?;
if reuse_would_collide(&chain, rest_ops) {
return None;
}
if reuse_would_create_back_edge(rank_of, &chain, rest_ops) {
return None;
}
Some(commit_reuse_chain(
nodes, edge_reads, fork, &chain, read_idx,
))
}
fn record_bypass_edge(
bypass_edges: &mut HashMap<usize, Vec<(usize, i32)>>,
from: usize,
to: usize,
) {
if from == to {
debug_assert_ne!(from, to, "record_bypass_edge: refusing a self-loop bypass");
return;
}
let entry = bypass_edges.entry(from).or_default();
for (t, w) in entry.iter_mut() {
if *t == to {
*w += 1;
return;
}
}
entry.push((to, 1));
}
#[allow(clippy::too_many_arguments)]
fn add_to_graph(
nodes: &mut Vec<Node>,
edge_reads: &mut HashMap<(usize, usize), Vec<u32>>,
bypass_edges: &mut HashMap<usize, Vec<(usize, i32)>>,
fork_arm_index: &mut HashMap<usize, HashMap<Vec<u8>, usize>>,
rank_of: &[usize],
query: &[u8],
ops: &[AlignOp],
read_idx: u32,
) {
let mut prev: Option<usize> = None;
let mut q_idx: usize = 0;
let mut i = 0usize;
let mut bypass_pending: Option<Option<usize>> = None;
let mut newly_forked: Vec<usize> = Vec::new();
while i < ops.len() {
let mut resume_node: Option<usize> = None;
let mut resume_is_bypass = false;
match ops[i] {
AlignOp::Match(node_idx) => {
let q_base = query[q_idx];
q_idx += 1;
let cur = if nodes[node_idx].base == q_base {
nodes[node_idx].coverage += 1;
if bypass_pending.is_some() {
resume_is_bypass = true;
} else if let Some(p) = prev {
if increment_or_add_edge(nodes, p, node_idx, false) {
newly_forked.push(p);
}
edge_reads.entry((p, node_idx)).or_default().push(read_idx);
}
node_idx
} else {
let edit = [q_base];
let reused = prev.and_then(|p| {
try_reuse_arm(
nodes,
edge_reads,
fork_arm_index,
rank_of,
p,
&edit,
&ops[i + 1..],
read_idx,
)
});
if let Some(reused_idx) = reused {
reused_idx
} else {
let new_idx = push_node(nodes, q_base);
nodes[new_idx].coverage = 1;
if let Some(p) = prev {
nodes[p].out_edges.push((
new_idx,
EdgeWeight {
matched: 1,
deleted: 0,
},
));
nodes[new_idx].in_edges.push(p);
edge_reads.entry((p, new_idx)).or_default().push(read_idx);
fork_arm_index
.entry(p)
.or_default()
.insert(edit.to_vec(), new_idx);
set_new_node_own_fork(nodes, p, new_idx);
newly_forked.push(p);
}
new_idx
}
};
resume_node = Some(cur);
prev = Some(cur);
i += 1;
}
AlignOp::Insert(first_base) => {
let run_start = i;
let mut edit = vec![first_base];
let mut j = i + 1;
while j < ops.len() {
if let AlignOp::Insert(b) = ops[j] {
edit.push(b);
j += 1;
} else {
break;
}
}
let reused = prev.and_then(|p| {
try_reuse_arm(
nodes,
edge_reads,
fork_arm_index,
rank_of,
p,
&edit,
&ops[j..],
read_idx,
)
});
if let Some(reused_idx) = reused {
resume_node = Some(reused_idx);
prev = Some(reused_idx);
} else {
let fork = prev;
let mut chain_start = None;
for &b in &edit {
let new_idx = push_node(nodes, b);
nodes[new_idx].coverage = 1;
if let Some(p) = prev {
nodes[p].out_edges.push((
new_idx,
EdgeWeight {
matched: 1,
deleted: 0,
},
));
nodes[new_idx].in_edges.push(p);
edge_reads.entry((p, new_idx)).or_default().push(read_idx);
set_new_node_own_fork(nodes, p, new_idx);
newly_forked.push(p);
}
chain_start.get_or_insert(new_idx);
prev = Some(new_idx);
}
if let (Some(f), Some(start)) = (fork, chain_start) {
fork_arm_index
.entry(f)
.or_default()
.insert(edit.clone(), start);
}
resume_node = chain_start;
}
q_idx += edit.len();
i = run_start + edit.len();
let _ = j; }
AlignOp::Delete(node_idx) => {
if bypass_pending.is_none() {
bypass_pending = Some(prev);
}
nodes[node_idx].delete_count += 1;
i += 1;
}
}
if bypass_pending.is_some() && resume_node.is_some() {
if resume_is_bypass {
if let (Some(Some(from)), Some(to)) = (bypass_pending, resume_node) {
debug_assert!(
from >= rank_of.len() || to >= rank_of.len() || rank_of[from] < rank_of[to],
"bypass edge {from}->{to} between in-snapshot nodes must respect \
topological order -- it is redundant with the real \
entry-pred->...->resume path through the skipped nodes"
);
debug_assert!(
edge_reads
.get(&(from, to))
.is_none_or(|v| v.last() != Some(&read_idx)),
"laundering guard: deleting read {read_idx} created a matched \
edge {from}->{to} into its bypass resume node"
);
record_bypass_edge(bypass_edges, from, to);
}
}
bypass_pending = None;
}
}
let _ = bypass_pending;
newly_forked.sort_unstable();
newly_forked.dedup();
for p in newly_forked {
propagate_fork_if_new(nodes, p);
}
}
fn heaviest_path(
nodes: &[Node],
topo: &[usize],
rank_of: &[usize],
bypass_edges: &HashMap<usize, Vec<(usize, i32)>>,
) -> Vec<(usize, u8, i32)> {
let n = topo.len();
let mut cum: Vec<(i64, Option<usize>, i32)> = vec![(0, None, 0); n];
for t in 0..n {
let node_idx = topo[t];
let node = &nodes[node_idx];
let curr = cum[t].0;
for &(succ_idx, ew) in &node.out_edges {
let succ_t = rank_of[succ_idx];
let candidate = curr + (ew.matched - 1) as i64;
if candidate > cum[succ_t].0 {
cum[succ_t] = (candidate, Some(t), ew.matched);
}
}
if let Some(bypasses) = bypass_edges.get(&node_idx) {
for &(succ_idx, weight) in bypasses {
let succ_t = rank_of[succ_idx];
debug_assert!(
t < succ_t,
"bypass edge {node_idx}->{succ_idx} (rank {t}->{succ_t}) must respect \
topological order; it is redundant with the real \
from->...->to path through the skipped nodes, so a violation means \
the graph is not a DAG"
);
let candidate = curr + (weight - 1) as i64;
if candidate > cum[succ_t].0 {
cum[succ_t] = (candidate, Some(t), weight);
}
}
}
}
let max_cum = (0..n).map(|t| cum[t].0).max().unwrap_or(0);
let best_t = (0..n).find(|&t| cum[t].0 == max_cum).unwrap_or(0);
let mut path: Vec<(usize, u8, i32)> = Vec::new();
let mut t = best_t;
loop {
let node_idx = topo[t];
let w = if cum[t].1.is_none() {
nodes[node_idx].coverage as i32
} else {
cum[t].2
};
path.push((node_idx, nodes[node_idx].base, w));
match cum[t].1 {
None => break,
Some(pred_t) => t = pred_t,
}
}
path.reverse();
path
}
fn greedy_heaviest_walk(
nodes: &[Node],
topo: &[usize],
rank_of: &[usize],
bypass_edges: &HashMap<usize, Vec<(usize, i32)>>,
) -> Vec<(usize, u8, i32)> {
let n = topo.len();
let mut score = vec![0i64; n];
let mut max_child = vec![usize::MAX; n];
let mut chosen_w = vec![0i32; n];
for t in (0..n).rev() {
let node_idx = topo[t];
let out_candidates = nodes[node_idx]
.out_edges
.iter()
.map(|&(succ, ew)| (rank_of[succ], ew.matched));
let bypass_candidates = bypass_edges
.get(&node_idx)
.into_iter()
.flatten()
.map(|&(succ, w)| (rank_of[succ], w));
let mut best: Option<(i32, usize)> = None; for (child_t, w) in out_candidates.chain(bypass_candidates) {
let better = match best {
None => true,
Some((bw, bc)) => w > bw || (w == bw && score[child_t] > score[bc]),
};
if better {
best = Some((w, child_t));
}
}
if let Some((w, child_t)) = best {
score[t] = w as i64 + score[child_t];
max_child[t] = child_t;
chosen_w[t] = w;
}
}
let src_t = (0..n)
.find(|&t| nodes[topo[t]].in_edges.is_empty())
.unwrap_or(0);
let mut path = Vec::new();
let mut t = src_t;
let mut incoming_w = nodes[topo[src_t]].coverage as i32;
loop {
let node_idx = topo[t];
path.push((node_idx, nodes[node_idx].base, incoming_w));
if max_child[t] == usize::MAX {
break;
}
incoming_w = chosen_w[t];
t = max_child[t];
}
path
}
fn majority_frequency(nodes: &[Node], topo: &[usize], min_cov: u32) -> Vec<(usize, u8, i32)> {
topo.iter()
.copied()
.filter(|&idx| {
let cov = nodes[idx].coverage;
let del = nodes[idx].delete_count;
let total = cov + del;
total >= min_cov && cov * 2 >= total
})
.map(|idx| (idx, nodes[idx].base, nodes[idx].coverage as i32))
.collect()
}
fn compute_stats(nodes: &[Node], min_allele_freq: f64, n_reads: usize) -> GraphStats {
let node_count = nodes.len();
let edge_count: usize = nodes.iter().map(|nd| nd.out_edges.len()).sum();
let coverages: Vec<f64> = nodes.iter().map(|nd| nd.coverage as f64).collect();
let coverage_mean = if node_count == 0 {
0.0
} else {
coverages.iter().sum::<f64>() / node_count as f64
};
let coverage_variance = if node_count == 0 {
0.0
} else {
coverages
.iter()
.map(|&c| (c - coverage_mean).powi(2))
.sum::<f64>()
/ node_count as f64
};
let single_support = nodes.iter().filter(|nd| nd.coverage == 1).count();
let single_support_fraction = if node_count == 0 {
0.0
} else {
single_support as f64 / node_count as f64
};
let mut weights: Vec<f64> = nodes
.iter()
.flat_map(|nd| nd.out_edges.iter().map(|&(_, ew)| ew.total() as f64))
.collect();
weights.sort_by(|a, b| a.partial_cmp(b).unwrap());
let edge_weight_gini = if weights.len() < 2 {
0.0
} else {
let n = weights.len() as f64;
let sum: f64 = weights.iter().sum();
if sum == 0.0 {
0.0
} else {
let numerator: f64 = weights
.iter()
.enumerate()
.map(|(i, &w)| (2.0 * (i as f64 + 1.0) - n - 1.0) * w)
.sum::<f64>();
numerator / (n * sum)
}
};
let threshold = (n_reads as f64 * min_allele_freq).ceil() as i32;
let mut bubble_count = 0usize;
let mut max_bubble_depth = 0usize;
let mut longest_bubble_span = 0usize;
for nd in nodes {
let qualifying: Vec<(usize, i32)> = nd
.out_edges
.iter()
.filter(|&&(_, ew)| ew.total() >= threshold)
.map(|&(to, ew)| (to, ew.total()))
.collect();
if qualifying.len() >= 2 {
bubble_count += 1;
let mut weights: Vec<i32> = qualifying.iter().map(|&(_, w)| w).collect();
weights.sort_unstable_by(|a, b| b.cmp(a));
max_bubble_depth = max_bubble_depth.max(weights[1] as usize);
for &(arm_start, _) in &qualifying {
let span = if nodes[arm_start].in_edges.len() > 1 {
0 } else {
materialize_arm_len(nodes, arm_start, ARM_MAX_DEPTH)
};
longest_bubble_span = longest_bubble_span.max(span);
}
}
}
let mean_column_entropy = {
let mut sum = 0.0f64;
let mut count = 0usize;
for nd in nodes {
let cov = nd.coverage as f64;
let del = nd.delete_count as f64;
let total = cov + del;
if total > 0.0 {
let h = binary_entropy(cov / total);
sum += h;
count += 1;
}
}
if count == 0 { 0.0 } else { sum / count as f64 }
};
GraphStats {
node_count,
edge_count,
bubble_count,
max_bubble_depth,
coverage_mean,
coverage_variance,
edge_weight_gini,
single_support_fraction,
mean_column_entropy,
longest_bubble_span,
median_input_read_len: 0, }
}
fn median_read_len(reads: &[Vec<u8>]) -> usize {
if reads.is_empty() {
return 0;
}
let mut lens: Vec<usize> = reads.iter().map(|r| r.len()).collect();
lens.sort_unstable();
lens[lens.len() / 2]
}
#[inline]
fn binary_entropy(p: f64) -> f64 {
if p <= 0.0 || p >= 1.0 {
0.0
} else {
let q = 1.0 - p;
-(p * p.log2() + q * q.log2())
}
}
fn detect_coverage_gaps(coverage: &[u32]) -> Vec<CoverageGap> {
let first = coverage.iter().position(|&c| c >= 2);
let last = coverage.iter().rposition(|&c| c >= 2);
let (first, last) = match (first, last) {
(Some(f), Some(l)) if f < l => (f, l),
_ => return vec![],
};
let mut gaps = Vec::new();
let mut gap_start: Option<usize> = None;
for (offset, &cov) in coverage[(first + 1)..last].iter().enumerate() {
let i = first + 1 + offset;
if cov < 2 {
gap_start.get_or_insert(i);
} else if let Some(s) = gap_start.take() {
gaps.push(CoverageGap {
start: s,
end: i,
kind: GapKind::Spanning,
});
}
}
if let Some(s) = gap_start {
gaps.push(CoverageGap {
start: s,
end: last,
kind: GapKind::Spanning,
});
}
gaps
}
fn find_bubbles(
nodes: &[Node],
topo: &[usize],
n_reads: usize,
min_allele_freq: f64,
) -> Vec<(usize, Vec<usize>)> {
let threshold = ((n_reads as f64 * min_allele_freq).ceil() as i32).max(1);
topo.iter()
.copied()
.filter_map(|node_idx| {
let arms: Vec<usize> = nodes[node_idx]
.out_edges
.iter()
.filter(|&&(_, ew)| ew.matched >= threshold)
.map(|&(to, _)| to)
.collect();
if arms.len() >= 2 {
Some((node_idx, arms))
} else {
None
}
})
.collect()
}
fn find_structural_bubbles(
nodes: &[Node],
topo: &[usize],
n_reads: usize,
cfg: &PoaConfig,
) -> Vec<(usize, Vec<usize>)> {
let threshold = ((n_reads as f64 * cfg.min_allele_freq).ceil() as i32).max(1);
let max_check = cfg.phasing_bubble_min_span.saturating_add(1);
topo.iter()
.copied()
.filter_map(|entry_node| {
let arms: Vec<usize> = nodes[entry_node]
.out_edges
.iter()
.filter(|&&(_, ew)| ew.matched >= threshold)
.map(|&(to, _)| to)
.collect();
if arms.len() < 2 {
return None;
}
let max_span = arms
.iter()
.map(|&start| {
if real_in_edge_count(nodes, start, threshold) > 1 {
0
} else {
materialize_arm_len_tolerant(nodes, start, max_check, threshold)
}
})
.max()
.unwrap_or(0);
if max_span >= cfg.phasing_bubble_min_span {
Some((entry_node, arms))
} else {
None
}
})
.collect()
}
const ARM_SEQUENCE_CAP: usize = 256;
fn arm_sequence(nodes: &[Node], start: usize) -> Vec<u8> {
if nodes[start].in_edges.len() > 1 {
return vec![];
}
let mut seq = Vec::new();
let mut cur = start;
loop {
if seq.len() >= ARM_SEQUENCE_CAP {
return vec![];
}
seq.push(nodes[cur].base);
match nodes[cur].out_edges.as_slice() {
[(next, _)] => {
if nodes[*next].in_edges.len() > 1 {
break; }
cur = *next;
}
_ => break, }
}
seq
}
fn collect_bubble_sites(
nodes: &[Node],
topo: &[usize],
filtered: &[(usize, u8, i32)],
edge_reads: &HashMap<(usize, usize), Vec<u32>>,
min_allele_freq: f64,
n_reads: usize,
phasing_bubble_min_span: usize,
) -> Vec<BubbleSite> {
let path_pos: HashMap<usize, usize> = filtered
.iter()
.enumerate()
.map(|(i, &(node_idx, _, _))| (node_idx, i))
.collect();
let threshold = ((n_reads as f64 * min_allele_freq).ceil() as i32).max(1);
find_bubbles(nodes, topo, n_reads, min_allele_freq)
.into_iter()
.filter_map(|(entry_node, arm_starts)| {
let consensus_pos = *path_pos.get(&entry_node)?;
let arm_read_counts: Vec<u32> = arm_starts
.iter()
.map(|&a| {
edge_reads
.get(&(entry_node, a))
.map_or(0, |v| v.len() as u32)
})
.collect();
let is_structural = arm_starts.iter().any(|&a| {
let len = if real_in_edge_count(nodes, a, threshold) > 1 {
0
} else {
materialize_arm_len_tolerant(nodes, a, phasing_bubble_min_span, threshold)
};
len >= phasing_bubble_min_span
});
let arm_sequences: Vec<Vec<u8>> =
arm_starts.iter().map(|&a| arm_sequence(nodes, a)).collect();
Some(BubbleSite {
consensus_pos,
arm_read_counts,
arm_sequences,
is_structural,
})
})
.collect()
}
const MIN_LENGTH_GAP_BP: f64 = 8.0;
const LENGTH_SEPARATION_MADS: f64 = 3.0;
const MIN_SPREAD_FLOOR_BP: f64 = 3.0;
fn median_usize(vals: &[usize]) -> f64 {
if vals.is_empty() {
return 0.0;
}
let mut v: Vec<f64> = vals.iter().map(|&x| x as f64).collect();
v.sort_by(|a, b| a.partial_cmp(b).unwrap());
let n = v.len();
if n % 2 == 1 {
v[n / 2]
} else {
(v[n / 2 - 1] + v[n / 2]) / 2.0
}
}
fn mad_usize(vals: &[usize], median: f64) -> f64 {
let dev: Vec<usize> = vals
.iter()
.map(|&x| (x as f64 - median).abs().round() as usize)
.collect();
median_usize(&dev)
}
fn length_separated(lens_a: &[usize], lens_b: &[usize]) -> bool {
if lens_a.is_empty() || lens_b.is_empty() {
return false;
}
let med_a = median_usize(lens_a);
let med_b = median_usize(lens_b);
let gap = (med_a - med_b).abs();
if gap < MIN_LENGTH_GAP_BP {
return false;
}
let spread_a = mad_usize(lens_a, med_a).max(MIN_SPREAD_FLOOR_BP);
let spread_b = mad_usize(lens_b, med_b).max(MIN_SPREAD_FLOOR_BP);
gap >= LENGTH_SEPARATION_MADS * spread_a.max(spread_b)
}
fn validate_and_merge_groups(
groups: Vec<Vec<usize>>,
reads: &[Vec<u8>],
min_reads: usize,
bubbles: &[(usize, Vec<usize>)],
edge_reads: &HashMap<(usize, usize), Vec<u32>>,
) -> Vec<Vec<usize>> {
if groups.len() < 2 {
return groups;
}
let mut groups = groups;
groups.sort_unstable_by_key(|g| std::cmp::Reverse(g.len()));
let n_bubbles = bubbles.len();
let require_length_separation = n_bubbles >= 2;
const CLEAR_MAJORITY: f64 = 0.60;
const MIN_ARM_COV: usize = 2;
let arm_of: Vec<HashMap<usize, usize>> = bubbles
.iter()
.map(|(entry, arm_starts)| {
let mut m: HashMap<usize, usize> = HashMap::new();
for (k, &a) in arm_starts.iter().enumerate() {
if let Some(rs) = edge_reads.get(&(*entry, a)) {
for &r in rs {
m.entry(r as usize).or_insert(k);
}
}
}
m
})
.collect();
let group_bubble_stats = |grp: &[usize], b: usize| -> (Option<usize>, f64, usize) {
let mut counts: HashMap<usize, usize> = HashMap::new();
for &r in grp {
if let Some(&k) = arm_of[b].get(&r) {
*counts.entry(k).or_default() += 1;
}
}
let cov: usize = counts.values().sum();
match counts.into_iter().max_by_key(|&(_, c)| c) {
Some((arm, c)) => (Some(arm), c as f64 / cov as f64, cov),
None => (None, 0.0, 0),
}
};
let clean_distinguishing = |ga: &[usize], gb: &[usize]| -> bool {
(0..n_bubbles).any(|b| {
let (ma, fa, ca) = group_bubble_stats(ga, b);
let (mb, fb, cb) = group_bubble_stats(gb, b);
ma.is_some()
&& mb.is_some()
&& ma != mb
&& fa >= CLEAR_MAJORITY
&& fb >= CLEAR_MAJORITY
&& ca >= MIN_ARM_COV
&& cb >= MIN_ARM_COV
})
};
let mut confirmed: Vec<Vec<usize>> = vec![groups[0].clone()];
let mut confirmed_lens: Vec<Vec<usize>> =
vec![groups[0].iter().map(|&r| reads[r].len()).collect()];
for cand in groups.into_iter().skip(1) {
let cand_lens: Vec<usize> = cand.iter().map(|&r| reads[r].len()).collect();
let significant = cand.len() >= min_reads;
let separated_from_all = !require_length_separation
|| confirmed_lens
.iter()
.all(|existing| length_separated(&cand_lens, existing));
if significant && separated_from_all {
let structurally_distinct = !require_length_separation
|| confirmed.iter().all(|c| clean_distinguishing(&cand, c));
if structurally_distinct {
confirmed.push(cand);
confirmed_lens.push(cand_lens);
} else {
let cand_med = median_usize(&cand_lens);
let target = (0..confirmed.len())
.filter(|&c| !clean_distinguishing(&cand, &confirmed[c]))
.min_by(|&a, &b| {
let da = (median_usize(&confirmed_lens[a]) - cand_med).abs();
let db = (median_usize(&confirmed_lens[b]) - cand_med).abs();
da.partial_cmp(&db).unwrap()
})
.expect("structurally_distinct false => at least one indistinguishable group");
confirmed_lens[target].extend(cand_lens);
confirmed[target].extend(cand);
}
} else if confirmed.len() == 1 {
confirmed_lens[0].extend(cand_lens);
confirmed[0].extend(cand);
} else {
let snapshot_medians: Vec<f64> =
confirmed_lens.iter().map(|l| median_usize(l)).collect();
let mut per_read_targets: Vec<(usize, usize)> = Vec::with_capacity(cand.len());
for (&r, &len) in cand.iter().zip(cand_lens.iter()) {
let len = len as f64;
let target = snapshot_medians
.iter()
.enumerate()
.min_by(|&(_, &a), &(_, &b)| {
let da = (a - len).abs();
let db = (b - len).abs();
da.partial_cmp(&db).unwrap()
})
.map(|(i, _)| i)
.unwrap_or(0);
per_read_targets.push((r, target));
}
for (r, target) in per_read_targets {
confirmed_lens[target].push(reads[r].len());
confirmed[target].push(r);
}
}
}
confirmed
}
fn phasing_groups(
edge_reads: &HashMap<(usize, usize), Vec<u32>>,
bubbles: &[(usize, Vec<usize>)],
n_reads: usize,
min_reads: usize,
reads: &[Vec<u8>],
) -> Vec<Vec<usize>> {
let n_bubbles = bubbles.len();
let mut sig: Vec<Vec<Option<usize>>> = vec![vec![None; n_bubbles]; n_reads];
for (b, (entry, arm_starts)) in bubbles.iter().enumerate() {
for (arm_idx, &arm_start) in arm_starts.iter().enumerate() {
if let Some(reads) = edge_reads.get(&(*entry, arm_start)) {
for &r in reads {
let r = r as usize;
if r < n_reads {
sig[r][b] = Some(arm_idx);
}
}
}
}
}
let mut assigned: Vec<usize> = Vec::new();
let mut unassigned: Vec<usize> = Vec::new();
for (r, row) in sig.iter().enumerate() {
if row.iter().any(|s| s.is_some()) {
assigned.push(r);
} else {
unassigned.push(r);
}
}
let mut order = assigned.clone();
order.sort_by_key(|&r| {
let known = sig[r].iter().filter(|s| s.is_some()).count();
(std::cmp::Reverse(known), r)
});
let mut clusters: Vec<Vec<Option<usize>>> = Vec::new();
let mut cluster_members: Vec<Vec<usize>> = Vec::new();
let mut bridge_candidates: Vec<(usize, Vec<usize>)> = Vec::new();
for r in order {
let row = &sig[r];
let compatible_clusters: Vec<usize> = clusters
.iter()
.enumerate()
.filter(|(_, known)| {
(0..n_bubbles).all(|b| match (row[b], known[b]) {
(Some(a), Some(bv)) => a == bv,
_ => true,
})
})
.map(|(ci, _)| ci)
.collect();
match compatible_clusters.as_slice() {
[] => {
clusters.push(row.clone());
cluster_members.push(vec![r]);
}
[only] => {
cluster_members[*only].push(r);
for b in 0..n_bubbles {
if clusters[*only][b].is_none() {
clusters[*only][b] = row[b];
}
}
}
_ => bridge_candidates.push((r, compatible_clusters)),
}
}
const PLAUSIBLE_LEN_FRACTION: f64 = 0.85;
if !cluster_members.is_empty() {
let snapshot_lens: Vec<f64> = cluster_members
.iter()
.map(|members| {
median_usize(&members.iter().map(|&r| reads[r].len()).collect::<Vec<_>>())
})
.collect();
let snapshot_sizes: Vec<usize> = cluster_members.iter().map(|m| m.len()).collect();
let min_full_len = (0..cluster_members.len())
.filter(|&c| snapshot_sizes[c] >= min_reads)
.map(|c| snapshot_lens[c])
.fold(f64::INFINITY, f64::min);
let nearest_cluster_among = |read_idx: usize, candidates: &[usize]| -> usize {
let len = reads[read_idx].len() as f64;
candidates
.iter()
.copied()
.min_by(|&a, &b| {
let da = (snapshot_lens[a] - len).abs();
let db = (snapshot_lens[b] - len).abs();
da.partial_cmp(&db)
.unwrap()
.then_with(|| snapshot_sizes[b].cmp(&snapshot_sizes[a]))
.then_with(|| a.cmp(&b))
})
.unwrap()
};
let largest_overall = (0..cluster_members.len())
.max_by_key(|&c| snapshot_sizes[c])
.unwrap();
let all_clusters: Vec<usize> = (0..cluster_members.len()).collect();
let largest_compatible = |candidates: &[usize]| -> usize {
*candidates
.iter()
.max_by_key(|&&c| snapshot_sizes[c])
.unwrap()
};
let mut resolved: Vec<(usize, usize)> = Vec::new();
for (r, compat) in &bridge_candidates {
let len = reads[*r].len() as f64;
let target = if min_full_len.is_finite() && len >= PLAUSIBLE_LEN_FRACTION * min_full_len
{
nearest_cluster_among(*r, &all_clusters)
} else {
largest_compatible(compat)
};
resolved.push((*r, target));
}
for &r in &unassigned {
let len = reads[r].len() as f64;
let target = if min_full_len.is_finite() && len >= PLAUSIBLE_LEN_FRACTION * min_full_len
{
nearest_cluster_among(r, &all_clusters)
} else {
largest_overall
};
resolved.push((r, target));
}
for (r, target) in resolved {
cluster_members[target].push(r);
}
} else {
cluster_members.push(Vec::new());
cluster_members[0].extend(bridge_candidates.into_iter().map(|(r, _)| r));
cluster_members[0].extend(unassigned);
}
let mut groups: Vec<Vec<usize>> = cluster_members;
groups.sort_unstable_by_key(|g| std::cmp::Reverse(g.len()));
if groups.is_empty() {
groups.push(Vec::new());
}
let mut i = 1;
while i < groups.len() {
if groups[i].len() < min_reads {
let g = groups.remove(i);
groups[0].extend(g);
} else {
i += 1;
}
}
groups.retain(|g| !g.is_empty());
groups
}
fn partition_reads_by_bubble(
edge_reads: &HashMap<(usize, usize), Vec<u32>>,
entry_node: usize,
arm_starts: &[usize],
n_reads: usize,
) -> Vec<Vec<usize>> {
let mut arm_order: Vec<usize> = (0..arm_starts.len()).collect();
arm_order.sort_unstable_by(|&a, &b| {
let wa = edge_reads
.get(&(entry_node, arm_starts[a]))
.map_or(0, |v| v.len());
let wb = edge_reads
.get(&(entry_node, arm_starts[b]))
.map_or(0, |v| v.len());
wb.cmp(&wa)
});
let n_arms = arm_order.len().min(2);
let mut groups: Vec<Vec<usize>> = vec![Vec::new(); n_arms];
let mut assigned = vec![false; n_reads];
for (slot, &arm_idx) in arm_order[..n_arms].iter().enumerate() {
let arm_start = arm_starts[arm_idx];
if let Some(reads) = edge_reads.get(&(entry_node, arm_start)) {
for &r in reads {
let r = r as usize;
if r < n_reads && !assigned[r] {
groups[slot].push(r);
assigned[r] = true;
}
}
}
}
let largest = groups
.iter()
.enumerate()
.max_by_key(|(_, g)| g.len())
.map(|(i, _)| i)
.unwrap_or(0);
for (r, &done) in assigned.iter().enumerate() {
if !done {
groups[largest].push(r);
}
}
groups.retain(|g| !g.is_empty());
groups
}
fn choose_seed(group: &[usize], reads: &[Vec<u8>]) -> usize {
if group.len() == 1 {
return 0;
}
let mut lens: Vec<usize> = group.iter().map(|&i| reads[i].len()).collect();
lens.sort_unstable();
let median = lens[lens.len() / 2];
group
.iter()
.enumerate()
.min_by_key(|&(_, &i)| reads[i].len().abs_diff(median))
.map(|(slot, _)| slot)
.unwrap_or(0)
}
impl PoaGraph {
pub fn new(seed: &[u8], config: PoaConfig) -> Result<Self, PoaError> {
if seed.is_empty() {
return Err(PoaError::EmptyInput);
}
let mut nodes: Vec<Node> = Vec::with_capacity(seed.len());
for &base in seed {
let idx = push_node(&mut nodes, base);
nodes[idx].coverage = 1;
}
let n = nodes.len();
let mut edge_reads: HashMap<(usize, usize), Vec<u32>> = HashMap::new();
for i in 0..n.saturating_sub(1) {
add_edge(&mut nodes, i, i + 1);
edge_reads.insert((i, i + 1), vec![0u32]);
}
Ok(PoaGraph {
nodes,
config,
n_reads: 1,
reads: vec![seed.to_vec()],
edge_reads,
edge_delete_reads: HashMap::new(),
bypass_edges: HashMap::new(),
warnings: 0,
cached_spine: Vec::new(),
spine_updated_at: 0,
spine_interval: 1,
align_scratch: AlignScratch::new(),
spine_mers: HashMap::new(),
fork_arm_index: HashMap::new(),
used_band_retry: false,
})
}
pub fn add_read(&mut self, read: &[u8]) -> Result<(), PoaError> {
if read.is_empty() {
return Err(PoaError::EmptyInput);
}
let (topo, rank_of) = topological_order(&self.nodes);
debug_assert_eq!(
topo.len(),
self.nodes.len(),
"graph has a cycle before adding read {}: topo includes {} of {} nodes",
self.n_reads,
topo.len(),
self.nodes.len()
);
let reads_since_update = self.n_reads.saturating_sub(self.spine_updated_at);
if self.cached_spine.is_empty() || reads_since_update >= self.spine_interval {
let new_spine = heaviest_path(&self.nodes, &topo, &rank_of, &self.bypass_edges);
let diff = spine_diff(&self.cached_spine, &new_spine);
self.cached_spine = new_spine;
self.spine_updated_at = self.n_reads;
self.spine_mers = build_spine_mers(&self.cached_spine, MINI_K, MINI_W);
if diff <= SPINE_STABLE_THRESHOLD {
self.spine_interval = (self.spine_interval * 2).min(SPINE_MAX_INTERVAL);
} else {
self.spine_interval = 1;
}
}
let read_mers = compute_minimizers(read, MINI_K, MINI_W);
let raw_anchors =
build_anchor_chain(&read_mers, &self.spine_mers, &self.cached_spine, &rank_of);
let anchors: Vec<(usize, usize)> = if raw_anchors.len() >= MINI_MIN_CHAIN {
raw_anchors
} else {
vec![]
};
let (ops, retried) = align_with_retry(
&self.nodes,
&topo,
&rank_of,
&self.cached_spine,
read,
&self.config,
&mut self.align_scratch,
&anchors,
)?;
if retried {
self.used_band_retry = true;
}
let read_idx = self.n_reads as u32;
add_to_graph(
&mut self.nodes,
&mut self.edge_reads,
&mut self.bypass_edges,
&mut self.fork_arm_index,
&rank_of,
read,
&ops,
read_idx,
);
self.reads.push(read.to_vec());
self.n_reads += 1;
Ok(())
}
pub fn consensus(&self) -> Result<Consensus, PoaError> {
if self.n_reads < self.config.min_reads {
return Err(PoaError::InsufficientDepth {
got: self.n_reads,
min: self.config.min_reads,
});
}
if self.n_reads == 1 {
let (topo, _) = topological_order(&self.nodes);
let sequence: Vec<u8> = topo.iter().map(|&idx| self.nodes[idx].base).collect();
let coverage: Vec<u32> = topo.iter().map(|_| 1).collect();
let path_weights: Vec<i32> = topo.iter().map(|_| 1).collect();
let mut graph_stats =
compute_stats(&self.nodes, self.config.min_allele_freq, self.n_reads);
graph_stats.median_input_read_len = median_read_len(&self.reads);
return Ok(Consensus {
sequence,
coverage,
path_weights,
n_reads: 1,
graph_stats,
gaps: vec![],
bubble_sites: vec![],
read_indices: vec![],
});
}
let (topo, rank_of) = topological_order(&self.nodes);
let filtered: Vec<(usize, u8, i32)> = match self.config.consensus_mode {
ConsensusMode::HeaviestPath => {
let path = greedy_heaviest_walk(&self.nodes, &topo, &rank_of, &self.bypass_edges);
let radius = LOCAL_POP_RADIUS;
let path_coverage: Vec<u32> = path
.iter()
.map(|&(node_idx, _, _)| self.nodes[node_idx].coverage)
.collect();
let local_pop = local_population_profile(&path_coverage, radius);
let local_min_cov_by_pos: Vec<u32> = local_pop
.iter()
.map(|&pop| coverage_threshold(pop as usize, self.config.min_coverage_fraction))
.collect();
let meets_floor: Vec<bool> = path
.iter()
.enumerate()
.map(|(i, &(node_idx, _, _))| {
self.nodes[node_idx].coverage >= local_min_cov_by_pos[i]
})
.collect();
let start = meets_floor.iter().position(|&s| s).unwrap_or(0);
let end = meets_floor
.iter()
.rposition(|&s| s)
.map(|i| i + 1)
.unwrap_or(path.len());
let range = if start < end {
start..end
} else {
0..path.len()
};
let range_start = range.start;
let range_start_rank = rank_of[path[range_start].0];
range
.filter(|&i| {
let (node_idx, _, _) = path[i];
if meets_floor[i] {
return true;
}
let fork_info = self.nodes[node_idx]
.nearest_fork
.filter(|&(pred_idx, _)| rank_of[pred_idx] >= range_start_rank);
let Some((pred_idx, arm_idx)) = fork_info else {
return self.nodes[node_idx].coverage
> self.nodes[node_idx].delete_count;
};
let local_total: i32 = self.nodes[pred_idx]
.out_edges
.iter()
.map(|&(_, ew)| ew.total())
.sum();
if (local_total.max(0) as u32) < local_min_cov_by_pos[i] {
return false;
}
let local_min_cov = coverage_threshold(
local_total.max(0) as usize,
self.config.min_coverage_fraction,
);
if self.nodes[node_idx].coverage as i32 >= local_min_cov as i32 {
return true;
}
if self.nodes[pred_idx].delete_count != 0
|| self.nodes[node_idx].delete_count != 0
{
return false;
}
let arm_weight = self.nodes[pred_idx]
.out_edges
.iter()
.find(|&&(to, _)| to == arm_idx)
.map(|&(_, ew)| ew.total())
.unwrap_or(0);
self.nodes[pred_idx]
.out_edges
.iter()
.all(|&(_, ew)| ew.total() <= arm_weight)
})
.map(|i| path[i])
.collect()
}
ConsensusMode::MajorityFrequency => {
majority_frequency(&self.nodes, &topo, self.min_coverage())
}
};
let sequence: Vec<u8> = filtered.iter().map(|&(_, base, _)| base).collect();
let coverage: Vec<u32> = filtered
.iter()
.map(|&(node_idx, _, _)| self.nodes[node_idx].coverage)
.collect();
let path_weights: Vec<i32> = filtered.iter().map(|&(_, _, w)| w).collect();
let mut graph_stats = compute_stats(&self.nodes, self.config.min_allele_freq, self.n_reads);
graph_stats.median_input_read_len = median_read_len(&self.reads);
let gaps = detect_coverage_gaps(&coverage);
let bubble_sites = collect_bubble_sites(
&self.nodes,
&topo,
&filtered,
&self.edge_reads,
self.config.min_allele_freq,
self.n_reads,
self.config.phasing_bubble_min_span,
);
Ok(Consensus {
sequence,
coverage,
path_weights,
n_reads: self.n_reads,
graph_stats,
gaps,
bubble_sites,
read_indices: vec![],
})
}
pub fn stats(&self) -> GraphStats {
let mut s = compute_stats(&self.nodes, self.config.min_allele_freq, self.n_reads);
s.median_input_read_len = median_read_len(&self.reads);
s
}
pub fn warnings_emitted(&self) -> usize {
self.warnings
}
pub fn edge_weights(&self) -> Vec<i32> {
self.nodes
.iter()
.flat_map(|n| n.out_edges.iter().map(|&(_, ew)| ew.total()))
.collect()
}
pub fn node_coverages(&self) -> Vec<u32> {
let (topo, _) = topological_order(&self.nodes);
topo.iter().map(|&i| self.nodes[i].coverage).collect()
}
pub fn align_read_ops(
&self,
read: &[u8],
) -> Result<(Vec<AlignOp>, usize, Vec<usize>), PoaError> {
let (topo, rank_of) = topological_order(&self.nodes);
let spine = heaviest_path(&self.nodes, &topo, &rank_of, &self.bypass_edges);
let (ops, _retried) = align_with_retry(
&self.nodes,
&topo,
&rank_of,
&spine,
read,
&self.config,
&mut AlignScratch::new(),
&[],
)?;
Ok((ops, 0, rank_of))
}
pub fn align_read_ops_unbanded(
&self,
read: &[u8],
) -> Result<(Vec<AlignOp>, Vec<usize>), PoaError> {
let (topo, rank_of) = topological_order(&self.nodes);
let spine = heaviest_path(&self.nodes, &topo, &rank_of, &self.bypass_edges);
let (ops, _retried) = align_with_retry(
&self.nodes,
&topo,
&rank_of,
&spine,
read,
&self.config,
&mut AlignScratch::new(),
&[],
)?;
Ok((ops, rank_of))
}
pub fn node_count(&self) -> usize {
self.nodes.len()
}
pub(crate) fn used_band_retry(&self) -> bool {
self.used_band_retry
}
pub fn graph_topology(&self) -> crate::types::GraphTopology {
use crate::types::{GraphEdgeInfo, GraphNodeInfo, GraphTopology};
let (topo, rank_of) = topological_order(&self.nodes);
let spine = heaviest_path(&self.nodes, &topo, &rank_of, &self.bypass_edges);
let nodes: Vec<GraphNodeInfo> = topo
.iter()
.enumerate()
.map(|(rank, &node_idx)| {
let n = &self.nodes[node_idx];
GraphNodeInfo {
node_idx,
base: n.base,
coverage: n.coverage,
delete_count: n.delete_count,
topo_rank: rank,
}
})
.collect();
let edges: Vec<GraphEdgeInfo> = topo
.iter()
.enumerate()
.flat_map(|(from_rank, &node_idx)| {
self.nodes[node_idx]
.out_edges
.iter()
.map(move |&(succ_idx, ew)| (from_rank, succ_idx, ew.total()))
.collect::<Vec<_>>()
})
.map(|(from_rank, succ_idx, weight)| GraphEdgeInfo {
from_rank,
to_rank: rank_of[succ_idx],
weight,
})
.collect();
let spine_ranks: Vec<usize> = spine
.iter()
.map(|(node_idx, _, _)| rank_of[*node_idx])
.collect();
GraphTopology {
nodes,
edges,
spine_ranks,
}
}
pub fn node_out_edge_info(&self) -> Vec<(usize, i32, i32)> {
self.nodes
.iter()
.map(|n| {
let count = n.out_edges.len();
let max_w = n
.out_edges
.iter()
.map(|&(_, ew)| ew.total())
.max()
.unwrap_or(0);
let min_w = n
.out_edges
.iter()
.map(|&(_, ew)| ew.total())
.min()
.unwrap_or(0);
(count, max_w, min_w)
})
.collect()
}
pub fn bubble_arm_lengths(
&self,
weight_threshold: i32,
arm_len_threshold: usize,
) -> Vec<(usize, Vec<usize>)> {
let (topo, _) = topological_order(&self.nodes);
let mut result = Vec::new();
for (t, &node_idx) in topo.iter().enumerate() {
let qualifying: Vec<usize> = self.nodes[node_idx]
.out_edges
.iter()
.filter(|&&(_, ew)| ew.total() >= weight_threshold)
.map(|&(succ, _)| succ)
.collect();
if qualifying.len() >= 2 {
let arm_lens: Vec<usize> = qualifying
.iter()
.map(|&arm_start| materialize_arm_len(&self.nodes, arm_start, 500))
.collect();
let min_len = arm_lens.iter().copied().min().unwrap_or(0);
if min_len >= arm_len_threshold {
result.push((t, arm_lens));
}
}
}
result
}
pub fn consensus_multi(&self) -> Result<Vec<Consensus>, PoaError> {
if self.n_reads < self.config.min_reads {
return Err(PoaError::InsufficientDepth {
got: self.n_reads,
min: self.config.min_reads,
});
}
let (topo, _) = topological_order(&self.nodes);
let structural = find_structural_bubbles(&self.nodes, &topo, self.n_reads, &self.config);
if structural.is_empty() && (self.config.band_width > 0 || self.config.adaptive_band) {
let mut cfg2 = self.config.clone();
cfg2.band_width = 0;
cfg2.adaptive_band = false;
let seed = &self.reads[0];
let mut g2 = PoaGraph::new(seed, cfg2)?;
for read in self.reads.iter().skip(1) {
g2.add_read(read)?;
}
return g2.consensus_multi();
}
let groups = if !structural.is_empty() {
let g = phasing_groups(
&self.edge_reads,
&structural,
self.n_reads,
self.config.min_reads,
&self.reads,
);
validate_and_merge_groups(
g,
&self.reads,
self.config.min_reads,
&structural,
&self.edge_reads,
)
} else {
let snp_bubbles = find_bubbles(
&self.nodes,
&topo,
self.n_reads,
self.config.min_allele_freq,
);
if !snp_bubbles.is_empty() {
let (entry, arm_starts) = snp_bubbles
.iter()
.max_by_key(|(entry, arms)| {
arms.iter()
.filter_map(|&arm| self.edge_reads.get(&(*entry, arm)))
.map(|v| v.len())
.min()
.unwrap_or(0)
})
.unwrap();
partition_reads_by_bubble(&self.edge_reads, *entry, arm_starts, self.n_reads)
} else {
return Ok(vec![self.consensus()?]);
}
};
if groups.len() < 2 {
return Ok(vec![self.consensus()?]);
}
let mut results = Vec::with_capacity(groups.len());
for group in &groups {
if group.len() < self.config.min_reads {
return Err(PoaError::InsufficientDepth {
got: group.len(),
min: self.config.min_reads,
});
}
let seed_slot = choose_seed(group, &self.reads);
let seed = &self.reads[group[seed_slot]];
let mut sub_cfg = self.config.clone();
sub_cfg.multi_allele = false;
let mut sub = PoaGraph::new(seed, sub_cfg)?;
for (slot, &read_idx) in group.iter().enumerate() {
if slot == seed_slot {
continue;
}
sub.add_read(&self.reads[read_idx])?;
}
let mut c = sub.consensus()?;
c.read_indices = group.to_vec();
results.push(c);
}
Ok(results)
}
fn min_coverage(&self) -> u32 {
coverage_threshold(self.n_reads, self.config.min_coverage_fraction)
}
}
fn coverage_threshold(population: usize, min_coverage_fraction: f64) -> u32 {
if min_coverage_fraction > 0.0 {
((population as f64 * min_coverage_fraction).ceil() as u32).max(1)
} else if population <= 1 {
1
} else {
((population / 2 + 1).max(2)) as u32
}
}
fn local_population_profile(coverages: &[u32], radius: usize) -> Vec<u32> {
let n = coverages.len();
let mut out = vec![0u32; n];
if n == 0 {
return out;
}
let mut deque: std::collections::VecDeque<usize> = std::collections::VecDeque::new();
let mut right = 0usize;
for (i, out_i) in out.iter_mut().enumerate() {
let hi = (i + radius).min(n - 1);
while right <= hi {
while let Some(&back) = deque.back() {
if coverages[back] <= coverages[right] {
deque.pop_back();
} else {
break;
}
}
deque.push_back(right);
right += 1;
}
let lo = i.saturating_sub(radius);
while let Some(&front) = deque.front() {
if front < lo {
deque.pop_front();
} else {
break;
}
}
*out_i = deque.front().map(|&idx| coverages[idx]).unwrap_or(0);
}
out
}
#[cfg(test)]
mod fork_cache_tests {
use super::*;
fn b(s: &str) -> Vec<u8> {
s.as_bytes().to_vec()
}
#[test]
fn nearest_fork_updates_for_preexisting_descendants_after_late_fork() {
let seed = b("AAAACGTACGTACGTA");
let cfg = PoaConfig {
min_reads: 3,
band_width: 0,
adaptive_band: false,
warn_on_long_unbanded: false,
..Default::default()
};
let mut g = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..4 {
g.add_read(&seed).unwrap();
}
for idx in 5..seed.len() {
assert_eq!(
g.nodes[idx].nearest_fork, None,
"node {idx} should have no fork ancestor yet (pre-divergence)"
);
assert_eq!(
g.nodes[idx].in_edges.len(),
1,
"node {idx} should be a plain single-predecessor chain node"
);
}
let fork_node = 4; assert_eq!(g.nodes[fork_node].out_edges.len(), 1, "not yet a fork");
let orig_child = g.nodes[fork_node].out_edges[0].0;
assert_eq!(orig_child, 5);
let mut divergent = b("AAAAC");
divergent.push(b'T'); for _ in 0..3 {
g.add_read(&divergent).unwrap();
}
assert_eq!(
g.nodes[fork_node].out_edges.len(),
2,
"node {fork_node} should now be a genuine fork"
);
fn ground_truth_nearest_fork(nodes: &[Node], start: usize) -> Option<(usize, usize)> {
let mut cur = start;
loop {
if nodes[cur].in_edges.len() != 1 {
return None; }
let pred = nodes[cur].in_edges[0];
if nodes[pred].out_edges.len() >= 2 {
return Some((pred, cur));
}
cur = pred;
}
}
for idx in 5..seed.len() {
let expected = ground_truth_nearest_fork(&g.nodes, idx);
assert_eq!(
g.nodes[idx].nearest_fork, expected,
"node {idx}: cached nearest_fork should match the independently-computed \
ground truth after node {fork_node} became a fork post-hoc"
);
assert_eq!(
g.nodes[idx].nearest_fork,
Some((fork_node, orig_child)),
"node {idx}: expected the cache to have been updated by forward \
propagation to point at the newly-created fork"
);
}
}
#[test]
fn nearest_fork_same_read_reconvergence_is_not_stale_updated() {
let seed = b("AAAACGTACGTACGTA");
let cfg = PoaConfig {
min_reads: 3,
band_width: 0,
adaptive_band: false,
warn_on_long_unbanded: false,
..Default::default()
};
let mut g = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..4 {
g.add_read(&seed).unwrap();
}
let fork_node = 4;
let orig_child = g.nodes[fork_node].out_edges[0].0;
assert_eq!(orig_child, 5);
let mut divergent = b("AAAAC");
divergent.push(b'T'); divergent.extend_from_slice(&seed[6..]);
for _ in 0..3 {
g.add_read(&divergent).unwrap();
}
assert_eq!(
g.nodes[fork_node].out_edges.len(),
2,
"node {fork_node} should now be a genuine fork"
);
let reconv = 6usize;
assert_eq!(
g.nodes[reconv].in_edges.len(),
2,
"node {reconv} should have genuinely reconverged (2 predecessors) \
within the same read that created the fork"
);
fn ground_truth_nearest_fork(nodes: &[Node], start: usize) -> Option<(usize, usize)> {
let mut cur = start;
loop {
if nodes[cur].in_edges.len() != 1 {
return None;
}
let pred = nodes[cur].in_edges[0];
if nodes[pred].out_edges.len() >= 2 {
return Some((pred, cur));
}
cur = pred;
}
}
assert_eq!(
g.nodes[5].nearest_fork,
Some((fork_node, orig_child)),
"node 5 (between the fork and the reconvergence) should still be updated"
);
assert_eq!(
g.nodes[5].nearest_fork,
ground_truth_nearest_fork(&g.nodes, 5)
);
for idx in 6..seed.len() {
let expected = ground_truth_nearest_fork(&g.nodes, idx);
assert_eq!(
expected, None,
"sanity check on the ground-truth helper itself: node {idx} \
is past a real reconvergence, so ground truth must be None"
);
assert_eq!(
g.nodes[idx].nearest_fork, expected,
"node {idx}: must not be stale-updated across the reconvergence \
at node {reconv} -- this is the exact hazard deferred (end-of-read) \
propagation is meant to avoid"
);
}
}
}
#[cfg(test)]
mod bypass_edge_tests {
use super::*;
fn b(s: &str) -> Vec<u8> {
s.as_bytes().to_vec()
}
const SEED: &str = "ACGTGCATCGTAGCTA";
fn cfg_global_unbanded() -> PoaConfig {
PoaConfig {
band_width: 0,
adaptive_band: false,
alignment_mode: AlignmentMode::Global,
warn_on_long_unbanded: false,
min_reads: 2,
..Default::default()
}
}
fn flatten(g: &PoaGraph) -> Vec<(usize, usize, i32)> {
let mut out: Vec<(usize, usize, i32)> = g
.bypass_edges
.iter()
.flat_map(|(&from, tos)| tos.iter().map(move |&(to, w)| (from, to, w)))
.collect();
out.sort_unstable();
out
}
fn out_edge(g: &PoaGraph, from: usize, to: usize) -> Option<EdgeWeight> {
g.nodes[from]
.out_edges
.iter()
.find(|&&(t, _)| t == to)
.map(|&(_, ew)| ew)
}
fn assert_no_delete_bucket_or_delete_reads(g: &PoaGraph) {
for (idx, nd) in g.nodes.iter().enumerate() {
for &(to, ew) in &nd.out_edges {
assert_eq!(
ew.deleted, 0,
"pure bypass: out-edge {idx}->{to} must carry no deleted-bucket weight, \
got {}",
ew.deleted
);
}
}
assert!(
g.edge_delete_reads.is_empty(),
"pure bypass: edge_delete_reads must never be populated"
);
}
#[test]
fn single_base_deletion_records_one_bypass_edge() {
let mut g = PoaGraph::new(&b(SEED), cfg_global_unbanded()).unwrap();
let read = b("ACGTGCTCGTAGCTA"); g.add_read(&read).unwrap();
assert_eq!(
flatten(&g),
vec![(5, 7, 1)],
"one single-base deletion should record exactly one bypass edge 5->7 weight 1"
);
assert_eq!(
g.nodes[6].delete_count, 1,
"delete_count on the skipped node must still be incremented"
);
assert_eq!(
out_edge(&g, 6, 7).map(|ew| ew.matched),
Some(1),
"edge 6->7 must NOT have gained matched weight from the deleting read \
(that would be the laundering this rework removes)"
);
assert_no_delete_bucket_or_delete_reads(&g);
}
#[test]
fn multi_base_deletion_run_records_one_spanning_bypass_edge() {
let mut g = PoaGraph::new(&b(SEED), cfg_global_unbanded()).unwrap();
let read = b("ACGTGCCGTAGCTA"); g.add_read(&read).unwrap();
assert_eq!(
flatten(&g),
vec![(5, 8, 1)],
"a 2-base deletion run should record exactly one bypass edge 5->8 weight 1"
);
assert_eq!(g.nodes[6].delete_count, 1);
assert_eq!(g.nodes[7].delete_count, 1);
assert_eq!(out_edge(&g, 7, 8).map(|ew| ew.matched), Some(1));
assert_no_delete_bucket_or_delete_reads(&g);
}
#[test]
fn deletion_run_spanning_a_preexisting_fork() {
let mut cfg = cfg_global_unbanded();
cfg.multi_allele = true;
let mut g = PoaGraph::new(&b(SEED), cfg).unwrap();
let ins = b("ACGTGCATAACGTAGCTA"); g.add_read(&ins).unwrap();
assert_eq!(
g.nodes[7].out_edges.len(),
2,
"node 7 should now be a fork (matched arm to C@8, inserted arm to the new 'A')"
);
assert!(
flatten(&g).is_empty(),
"an insertion-only read must not record any bypass edge"
);
let del = b("ACGTGCGTAGCTA"); g.add_read(&del).unwrap();
g.add_read(&del).unwrap();
assert_eq!(
flatten(&g),
vec![(5, 9, 2)],
"a deletion run spanning a fork should still record one bypass edge 5->9, \
incremented to weight 2 across the two identical deletion reads"
);
assert_eq!(g.nodes[6].delete_count, 2);
assert_eq!(g.nodes[7].delete_count, 2);
assert_eq!(g.nodes[8].delete_count, 2);
assert_eq!(out_edge(&g, 8, 9).map(|ew| ew.matched), Some(2));
assert_no_delete_bucket_or_delete_reads(&g);
}
#[test]
fn leading_delete_run_records_no_bypass_edge() {
let mut g = PoaGraph::new(&b(SEED), cfg_global_unbanded()).unwrap();
let suffix = b("ATCGTAGCTA"); g.add_read(&suffix).unwrap();
assert!(
flatten(&g).is_empty(),
"a leading delete run (no predecessor to bypass from) must record no bypass edge"
);
assert_eq!(
g.nodes[0].delete_count, 1,
"leading deleted nodes are still counted by the existing bookkeeping"
);
}
#[test]
fn trailing_terminal_delete_run_records_no_bypass_edge() {
let mut g = PoaGraph::new(&b(SEED), cfg_global_unbanded()).unwrap();
let (_topo, rank_of) = topological_order(&g.nodes);
let ops = vec![
AlignOp::Match(0),
AlignOp::Match(1),
AlignOp::Delete(2),
AlignOp::Delete(3),
];
let query = b("AC"); add_to_graph(
&mut g.nodes,
&mut g.edge_reads,
&mut g.bypass_edges,
&mut g.fork_arm_index,
&rank_of,
&query,
&ops,
1,
);
assert!(
flatten(&g).is_empty(),
"a trailing terminal delete run must not create a bypass edge to nowhere"
);
assert_eq!(g.nodes[2].delete_count, 1);
assert_eq!(g.nodes[3].delete_count, 1);
assert_no_delete_bucket_or_delete_reads(&g);
}
#[test]
fn no_deletions_records_no_bypass_edges() {
let mut g = PoaGraph::new(&b(SEED), cfg_global_unbanded()).unwrap();
g.add_read(&b(SEED)).unwrap();
g.add_read(&b(SEED)).unwrap();
assert!(
flatten(&g).is_empty(),
"reads with no Delete ops must leave bypass_edges empty"
);
}
}
#[cfg(test)]
mod validate_and_merge_tests {
use super::*;
fn read_of_len(n: usize) -> Vec<u8> {
vec![b'A'; n]
}
#[test]
fn merges_structurally_indistinct_length_split() {
let mut reads = Vec::new();
for _ in 0..10 {
reads.push(read_of_len(100));
}
for k in 0..5 {
reads.push(read_of_len(314 + k)); }
for k in 0..4 {
reads.push(read_of_len(298 + k)); }
let g0: Vec<usize> = (0..10).collect();
let sub_a: Vec<usize> = (10..15).collect();
let sub_b: Vec<usize> = (15..19).collect();
let bubbles = vec![
(0usize, vec![1usize, 2usize]),
(3usize, vec![4usize, 5usize]),
];
let mut edge_reads: HashMap<(usize, usize), Vec<u32>> = HashMap::new();
edge_reads.insert((0, 1), g0.iter().map(|&r| r as u32).collect());
edge_reads.insert((3, 4), g0.iter().map(|&r| r as u32).collect());
let long: Vec<u32> = sub_a.iter().chain(&sub_b).map(|&r| r as u32).collect();
edge_reads.insert((0, 2), long.clone());
edge_reads.insert((3, 5), long);
let out = validate_and_merge_groups(
vec![g0.clone(), sub_a.clone(), sub_b.clone()],
&reads,
3,
&bubbles,
&edge_reads,
);
assert_eq!(
out.len(),
2,
"structurally-indistinct length-split sub-groups must merge to 2 alleles, got {}",
out.len()
);
let long_group = out
.iter()
.find(|g| g.contains(&10))
.expect("a group containing the long allele");
for r in sub_a.iter().chain(&sub_b) {
assert!(
long_group.contains(r),
"long-allele read {r} must be in the single merged long group"
);
}
}
#[test]
fn keeps_structurally_distinct_length_split() {
let mut reads = Vec::new();
for _ in 0..10 {
reads.push(read_of_len(100)); }
for _ in 0..8 {
reads.push(read_of_len(200)); }
let gx: Vec<usize> = (0..10).collect();
let gy: Vec<usize> = (10..18).collect();
let bubbles = vec![
(0usize, vec![1usize, 2usize]),
(3usize, vec![4usize, 5usize]),
];
let mut edge_reads: HashMap<(usize, usize), Vec<u32>> = HashMap::new();
edge_reads.insert((0, 1), gx.iter().map(|&r| r as u32).collect());
edge_reads.insert((3, 4), gx.iter().map(|&r| r as u32).collect());
edge_reads.insert((0, 2), gy.iter().map(|&r| r as u32).collect());
edge_reads.insert((3, 5), gy.iter().map(|&r| r as u32).collect());
let out = validate_and_merge_groups(vec![gx, gy], &reads, 3, &bubbles, &edge_reads);
assert_eq!(
out.len(),
2,
"two genuine alleles with a clean distinguishing bubble must stay split, got {}",
out.len()
);
}
}
#[cfg(test)]
mod band_too_narrow_tests {
use super::*;
fn fixture() -> (Vec<u8>, Vec<u8>) {
let seed: Vec<u8> = "ACGTACGTCG".repeat(10).into_bytes(); let mut query = seed.clone();
query.extend_from_slice(&"TGCA".repeat(30).into_bytes()); (seed, query)
}
fn align_direct(
g: &PoaGraph,
topo: &[usize],
rank_of: &[usize],
spine: &[(usize, u8, i32)],
query: &[u8],
band_width: usize,
scratch: &mut AlignScratch,
) -> Result<Vec<AlignOp>, PoaError> {
let cfg = PoaConfig {
band_width,
adaptive_band: false,
alignment_mode: AlignmentMode::SemiGlobal,
..PoaConfig::default()
};
align(&g.nodes, topo, rank_of, spine, query, &cfg, scratch, &[])
}
#[test]
fn band_too_narrow_returned_below_threshold_ok_at_and_above() {
let (seed, query) = fixture();
let g = PoaGraph::new(&seed, PoaConfig::default()).unwrap();
let (topo, rank_of) = topological_order(&g.nodes);
let spine = heaviest_path(&g.nodes, &topo, &rank_of, &g.bypass_edges);
let mut scratch = AlignScratch::new();
for &band_width in &[10usize, 50, 100] {
let result = align_direct(
&g,
&topo,
&rank_of,
&spine,
&query,
band_width,
&mut scratch,
);
match result {
Err(PoaError::BandTooNarrow {
configured,
required,
}) => {
assert_eq!(configured, band_width);
assert!(
required > band_width,
"required ({required}) should exceed the too-narrow configured \
width ({band_width})"
);
}
other => panic!(
"band_width={band_width} expected BandTooNarrow, got {:?}",
other.map(|ops| ops.len())
),
}
}
for &band_width in &[119usize, 120, 121, 150] {
let result = align_direct(
&g,
&topo,
&rank_of,
&spine,
&query,
band_width,
&mut scratch,
);
match result {
Ok(ops) => assert_eq!(
ops.len(),
query.len(),
"band_width={band_width}: expected one op per query base \
(100 Match + 120 Insert), got {} ops",
ops.len()
),
Err(e) => panic!("band_width={band_width} expected Ok, got {e:?}"),
}
}
}
#[test]
fn align_with_retry_recovers_correct_alignment_from_any_narrow_start() {
let (seed, query) = fixture();
let g = PoaGraph::new(&seed, PoaConfig::default()).unwrap();
let (topo, rank_of) = topological_order(&g.nodes);
let spine = heaviest_path(&g.nodes, &topo, &rank_of, &g.bypass_edges);
let mut scratch = AlignScratch::new();
for &band_width in &[10usize, 50, 100] {
let cfg = PoaConfig {
band_width,
adaptive_band: false,
alignment_mode: AlignmentMode::SemiGlobal,
..PoaConfig::default()
};
let (ops, retried) = align_with_retry(
&g.nodes,
&topo,
&rank_of,
&spine,
&query,
&cfg,
&mut scratch,
&[],
)
.unwrap_or_else(|e| panic!("band_width={band_width}: retry should recover, got {e:?}"));
assert!(
retried,
"band_width={band_width}: expected align_with_retry to report that pass 1 \
needed widening"
);
let n_match = ops
.iter()
.filter(|o| matches!(o, AlignOp::Match(_)))
.count();
let n_insert = ops
.iter()
.filter(|o| matches!(o, AlignOp::Insert(_)))
.count();
let n_delete = ops
.iter()
.filter(|o| matches!(o, AlignOp::Delete(_)))
.count();
assert_eq!(
(n_match, n_insert, n_delete),
(100, 120, 0),
"band_width={band_width}: expected the retry to recover the exact \
alignment shape (100 Match + 120 Insert, 0 Delete), got \
Match={n_match} Insert={n_insert} Delete={n_delete}"
);
}
}
}