poa_consensus/

lib.rs

//! # poa-consensus
//!
//! A pure-Rust banded Partial Order Alignment (POA) library for building a
//! consensus sequence from a set of reads.
//!
//! POA builds a directed acyclic graph (DAG) from the reads, aligns each
//! subsequent read into the graph using affine-gap dynamic programming, and
//! extracts a consensus by following the heaviest (most-supported) path through
//! the graph.  It handles length variation between reads naturally: inserts and
//! deletions create separate branches in the graph, and the heaviest-path
//! extraction resolves them by read support rather than by fixed rules.
//!
//! ## When to use this crate
//!
//! This crate is optimised for **short to medium reads** (50 bp to ~20 kb with
//! banded DP) where the graph stays small enough for in-memory DP — short
//! tandem repeat (STR) loci, amplicon consensus, per-locus nanopore or HiFi
//! read sets.
//!
//! The nearest alternative on crates.io is
//! [`poasta`](https://crates.io/crates/poasta) (Broad Institute, pure Rust,
//! gap-affine A\* alignment).  `poasta` excels at larger graphs such as
//! bacterial genes and HLA loci.  For STR reads where graphs are small and
//! throughput across many loci matters, a well-tuned banded DP is faster and
//! simpler.  Both crates are pure Rust; neither wraps a C library.
//!
//! ## Quick start
//!
//! ### Functional API (single call)
//!
//! ```rust
//! use poa_consensus::{consensus, consensus_multi, PoaConfig};
//!
//! let reads: Vec<&[u8]> = vec![
//!     b"CATCATCAT",
//!     b"CATCATCAT",
//!     b"CATCGTCAT",
//!     b"CATCATCAT",
//! ];
//!
//! // Single-allele consensus; seed_idx=0 seeds the graph with reads[0].
//! let result = consensus(&reads, 0, &PoaConfig::default())?;
//! println!("{}", String::from_utf8_lossy(&result.sequence));
//!
//! // Multi-allele: returns one Consensus per detected allele.
//! let alleles = consensus_multi(&reads, 0, &PoaConfig::default())?;
//! # Ok::<(), poa_consensus::PoaError>(())
//! ```
//!
//! ### Stateful API (inspect graph between reads)
//!
//! ```rust
//! use poa_consensus::{PoaGraph, PoaConfig};
//!
//! let reads: &[&[u8]] = &[b"CATCATCAT", b"CATCATCAT", b"CATCGTCAT"];
//!
//! let mut graph = PoaGraph::new(reads[0], PoaConfig::default())?;
//! for read in &reads[1..] {
//!     graph.add_read(read)?;
//! }
//! let consensus = graph.consensus()?;
//! let stats     = graph.stats();
//! println!("bubbles: {}", stats.bubble_count);
//! # Ok::<(), poa_consensus::PoaError>(())
//! ```
//!
//! ### Two-pass adaptive mode
//!
//! ```rust
//! use poa_consensus::{consensus_adaptive, PoaConfig};
//!
//! let reads: Vec<&[u8]> = vec![
//!     b"CATCATCAT", b"CATCATCAT", b"CATCATCAT",
//!     b"CATCGTCAT", b"CATCGTCAT", b"CATCGTCAT",
//! ];
//! // Pass 1 builds the graph and computes GraphStats.
//! // Pass 2 is selected automatically: multi-allele split, noise tightening,
//! // or semi-global switch, depending on what the stats reveal.
//! let result  = consensus_adaptive(&reads, 0, &PoaConfig::default())?;
//! let alleles = result.consensuses;   // Vec<Consensus>; one or two elements
//! # Ok::<(), poa_consensus::PoaError>(())
//! ```
//!
//! ## Seed selection
//!
//! The seed read initialises the graph as a linear chain of nodes.  All other
//! reads are aligned into this initial structure, so a poor seed degrades
//! alignment quality for every subsequent read.
//!
//! **Choose a median-length read.**  The seed acts as the backbone; reads
//! shorter than the seed produce terminal deletions and reads longer than the
//! seed produce extensions.  A median-length seed minimises both.  In
//! high-throughput use (many loci), `reads.iter().enumerate().min_by_key(|(_, r)| r.len().abs_diff(median_len))` is a one-liner.
//!
//! **Avoid outliers.**  The longest and shortest reads are the most likely to
//! be error-prone or to span a different number of repeat units.  Using one as
//! the seed skews the initial graph in a direction that the alignment of
//! subsequent reads must then correct.
//!
//! **Orient reads before selecting the seed.**  Mixed-strand input produces a
//! garbage graph silently.  Call [`auto_orient`] before POA construction; it
//! uses k-mer matching (O(n) per read, no alignment) and returns borrowed
//! slices for reads already on the correct strand.
//!
//! Seed selection is the caller's responsibility.  The API takes an explicit
//! `seed_idx` so that the selection logic can live in the caller and be tuned
//! per application.
//!
//! ## Band width and scale
//!
//! The DP matrix has one row per read base and one column per graph node.
//! Unbanded alignment is O(read_len × graph_nodes) memory and time — manageable
//! for reads up to ~1 kb but prohibitive above that.
//!
//! | Read length | Band width | Memory per read (3 matrices, i32) |
//! |---|---|---|
//! | 600 bp | unbanded | ~1.4 MB |
//! | 600 bp | 100 | ~2.4 KB |
//! | 20 kb | unbanded | ~9.6 GB |
//! | 20 kb | adaptive (w≈210) | ~200 MB |
//!
//! **Rules of thumb:**
//!
//! - Reads ≤ 1 kb: unbanded (`band_width = 0`) is fine.
//! - Reads 1 kb–20 kb: set `band_width` to at least twice the expected
//!   length difference between reads, or enable `adaptive_band = true` (abPOA
//!   formula: `w = adaptive_band_b + adaptive_band_f × read_len`).
//! - Reads > 20 kb: adaptive banding is required; consider `poasta` for graphs
//!   that approach bacterial-gene size.
//!
//! A band that is too narrow returns `Err(PoaError::BandTooNarrow)` when the
//! terminal column is unreachable.  The library never silently produces a wrong
//! alignment — it errors instead.  [`PoaGraph::warnings_emitted`] counts how
//! many times a long-unbanded warning fired during `add_read` calls.
//!
//! ## Coverage and depth
//!
//! **`min_reads`** (default: 1) is the minimum number of reads required to
//! attempt consensus.  `consensus()` returns `Err(InsufficientDepth)` below
//! this threshold.  For reliable results, use at least 5 reads; 10+ is
//! preferable for heterozygous sites.
//!
//! **Boundary trim** removes leading and trailing nodes whose coverage falls
//! below the majority threshold `(n/2 + 1).max(2)` (or
//! `min_coverage_fraction × n` if set explicitly).  This corrects for seed
//! reads that happen to be longer or shorter than the true consensus length.
//!
//! **In multi-allele mode**, `min_reads` is applied per allele group, not to
//! the total read count.  With 10 reads split 5/5, each group must independently
//! satisfy `min_reads`.
//!
//! ## Known limitations
//!
//! **Phase-shift majority trim.** When the majority of reads are phase-shifted
//! relative to the seed (e.g. most reads start one repeat unit later), boundary
//! trim may incorrectly remove the seed's first node because its Match coverage
//! is low.  The long-term fix is to use [`extract_flanked_region`] to anchor
//! reads to a common reference point before POA; this eliminates phase ambiguity
//! entirely.
//!
//! **Rotation-ambiguous repeats without flanking sequence.** Trinucleotide
//! repeats like GAA can appear as GAA, AAG, or AGA depending on where the read
//! starts relative to the period.  Without flanking anchors, POA on a mixed
//! rotational-phase read set produces unreliable output.  The fix is the same:
//! run [`extract_flanked_region`] first so that every read enters POA already
//! anchored at the same phase.
//!
//! **Long expansions and partial reads.** When the repeat is longer than most
//! reads, reads that do not span the full locus create artifactual terminal
//! deletions that do not increment node coverage, making boundary nodes appear
//! low-coverage and vulnerable to trim.  Use `AlignmentMode::SemiGlobal` (free
//! terminal gaps) to prevent partial reads from distorting the boundary, or
//! apply [`extract_flanked_region`] to exclude non-spanning reads before POA.
//!
//! ---
//!
//! ## Affine gap penalties
//!
//! The alignment scoring uses **affine gap penalties**. This section explains
//! what that means and why it matters, because the standard explanation in most
//! papers is needlessly opaque.
//!
//! ### The problem with linear gaps
//!
//! A linear gap model charges a flat penalty per gap base, say `-1` for every
//! missing or inserted base. This means a single deletion of 4 bases and four
//! separate single-base deletions cost exactly the same:
//!
//! ```text
//! One 4-base deletion:       -1 -1 -1 -1   = -4
//! Four 1-base deletions:     -1 -1 -1 -1   = -4
//! ```
//!
//! In biology these are very different events. A 4-base deletion is one
//! mutation: a single copying mistake that dropped 4 bases together. Four
//! separate 1-base deletions are four independent mutations. The linear model
//! treats them identically, so it cannot prefer the more plausible explanation.
//!
//! ### Affine gaps: one price to start, a smaller price to continue
//!
//! Affine scoring uses two parameters:
//!
//! - **`gap_open`**: a one-time penalty paid at the moment a gap starts.
//!   This represents the cost of the mutation event itself.
//! - **`gap_extend`**: a per-base penalty paid for every base inside the gap.
//!   This represents the cost of each missing or inserted base.
//!
//! A gap of length `k` costs: `gap_open + k * gap_extend`
//!
//! ### Example
//!
//! Using `gap_open = -2` and `gap_extend = -1`:
//!
//! ```text
//! One gap of length 4:
//!     gap_open + 4 * gap_extend = -2 + (-4) = -6
//!
//! Four gaps of length 1:
//!     4 * (gap_open + 1 * gap_extend) = 4 * (-3) = -12
//! ```
//!
//! The single 4-base deletion costs -6. Four separate 1-base deletions cost
//! -12. The aligner now strongly prefers the single-event explanation, which
//! matches biological reality.
//!
//! With linear gaps (gap = -1) both scenarios cost -4 and the aligner cannot
//! distinguish them at all.
//!
//! ### Worked alignment example
//!
//! Say we are aligning a read to a reference with a 4-base deletion:
//!
//! ```text
//! Reference: A C G T T T T T A C G T
//! Read:      A C G T - - - - A C G T
//! ```
//!
//! Scoring (match = +1, mismatch = -1, gap_open = -2, gap_extend = -1):
//!
//! ```text
//! 8 matching bases:  8 * (+1)           =  +8
//! 1 gap of length 4: -2 + 4 * (-1)     =  -6
//! Total:                                =  +2
//! ```
//!
//! Now consider a different alignment that avoids the gap by accepting 4
//! mismatches instead:
//!
//! ```text
//! Reference: A C G T T T T T A C G T
//! Read:      A C G T X X X X A C G T   (X = mismatch)
//! ```
//!
//! ```text
//! 8 matching bases:   8 * (+1)          =  +8
//! 4 mismatches:       4 * (-1)          =  -4
//! Total:                                =  +4
//! ```
//!
//! Here the mismatched alignment scores higher (+4 vs +2), so the aligner
//! would choose it. Whether that is the right call depends on the biology: if
//! the read truly has a 4-base deletion, these parameters would produce the
//! wrong alignment. Tuning `gap_open` and `gap_extend` controls this
//! trade-off. A more negative `gap_open` makes the aligner more willing to
//! open a gap rather than accept mismatches; a less negative `gap_extend`
//! makes longer gaps cheaper relative to short ones.
//!
//! ### Why individual alignment ambiguity matters less in POA
//!
//! In a single pairwise alignment there is often no way to know whether a gap
//! or a mismatch is the correct interpretation. POA largely sidesteps this
//! problem because it aligns many reads into the same graph and then extracts
//! a consensus, rather than making a definitive call on any one alignment in
//! isolation.
//!
//! If a gap in one read reflects a real deletion, most other reads will carry
//! the same gap. The corresponding graph edges accumulate high weight, and the
//! heaviest-path consensus will include the deletion. If the gap is a
//! sequencing error, other reads will traverse the same graph node as a match,
//! and the consensus will correct it.
//!
//! This means parameter choices affect alignment quality at the margins, but
//! the coverage threshold and heaviest-path extraction together resolve most
//! of the ambiguity that would be fatal to a single pairwise alignment. The
//! defaults (`gap_open = -2`, `gap_extend = -1`) are calibrated for Oxford
//! Nanopore and PacBio HiFi reads, where indels are the dominant error type.
//!
//! ### How affine scoring changes the dynamic programming
//!
//! Standard Needleman-Wunsch uses a single DP table. Affine scoring requires
//! three tables, one per "state" the aligner can be in at any position:
//!
//! - **`M[i][j]`**: best score where query position `i` aligns to graph node
//!   `j` as a match or mismatch.
//! - **`I[i][j]`**: best score where query position `i` is an insertion (a
//!   base in the read with no corresponding node in the graph).
//! - **`D[i][j]`**: best score where graph node `j` is a deletion (a node in
//!   the graph skipped by the read).
//!
//! The recurrences are:
//!
//! ```text
//! M[i][j] = max(M[i-1][j-1], I[i-1][j-1], D[i-1][j-1]) + score(query[i], node[j])
//!
//! I[i][j] = max(M[i][j-1] + gap_open + gap_extend,
//!               I[i][j-1] + gap_extend)
//!
//! D[i][j] = max(M[i-1][j] + gap_open + gap_extend,
//!               D[i-1][j] + gap_extend)
//! ```
//!
//! The key insight in `I` and `D`: transitioning from `M` to a gap state pays
//! `gap_open + gap_extend` (opening cost plus first base). Staying in the same
//! gap state pays only `gap_extend` (continuing an existing gap). The aligner
//! tracks which state it is in at every cell, so it knows whether to apply the
//! full open penalty or just the extend penalty.
//!
//! Traceback follows whichever state and predecessor gave the maximum score at
//! each cell, producing a sequence of Match, Insert, and Delete operations.
//!
//! For POA over a graph rather than a linear reference, `j` indexes a graph
//! node in topological order rather than a column in a matrix, but the
//! recurrence is identical.
//!
//! ### Choosing parameter values
//!
//! The defaults used in this crate (`gap_open = -2`, `gap_extend = -1`) work
//! well for Oxford Nanopore and PacBio HiFi reads at short tandem repeat loci.
//! As a rule of thumb:
//!
//! - Increase the magnitude of `gap_open` (e.g. `-4`) to make the aligner
//!   more reluctant to open gaps, preferring mismatches instead.
//! - Decrease the magnitude of `gap_extend` (e.g. `-0.5`, or use integer
//!   scaling) to make long gaps cheaper, useful when reads have systematic
//!   length variation.
//! - Setting `gap_open = 0` recovers linear gap behaviour where only
//!   `gap_extend` matters.

pub mod analysis;
pub mod config;
pub mod error;
pub mod flank;
pub mod graph;
pub mod orient;
pub mod seed;
pub mod types;

#[cfg(feature = "plot")]
pub mod plot;

#[cfg(test)]
mod tests;

pub use analysis::{
    ConsensusWarnings, DiagnoseConfig, InteriorSupportWarning, LowDepthWarning,
    StructuralCompetingSummary, TruncationWarning, diagnose,
};
pub use config::{AlignmentMode, ConsensusMode, PoaConfig};
pub use error::PoaError;
pub use flank::extract_flanked_region;
pub use graph::{AlignOp, PoaGraph};
pub use orient::{auto_orient, orient_to_seed, reverse_complement};
pub use seed::{SeedSelection, select_seed};
pub use types::{
    AdaptiveAction, AdaptiveResult, BubbleSite, Consensus, CoverageGap, GapKind, GraphEdgeInfo,
    GraphNodeInfo, GraphStats, GraphTopology, Strand,
};

// ── Internal helpers ─────────────────────────────────────────────────────────

fn build_graph(reads: &[&[u8]], seed_idx: usize, config: PoaConfig) -> Result<PoaGraph, PoaError> {
    let mut graph = PoaGraph::new(reads[seed_idx], config.clone())?;
    for (i, read) in reads.iter().enumerate() {
        if i != seed_idx {
            graph.add_read(read)?;
        }
    }

    // If any read needed a wider-than-configured band during the
    // incremental build (`align_with_retry`'s pass 2 or 3 -- see
    // `PoaGraph::used_band_retry`'s doc comment), the resulting graph is a
    // mix of different reads settling with different effective band
    // widths. In a periodic/repetitive locus this can settle different
    // reads onto different, individually-plausible diagonals (the same
    // mechanism as Known Bug #3's "silent wrong alignment on narrow band
    // in repetitive regions"), fragmenting bubble structure that a
    // uniformly-unbanded build would not -- confirmed on real data
    // (`multi_gaa30_100`): a mixed-band graph produced a spurious 3rd
    // allele that vanished when every read used the same unbanded config
    // from the start. Rebuild the whole graph unbanded, from scratch, for
    // internal consistency, rather than trusting a mixed-band graph as-is.
    // Skipped when the config was already fully unbanded (nothing to gain,
    // and this recursion would otherwise never terminate).
    if graph.used_band_retry() && (config.band_width > 0 || config.adaptive_band) {
        let mut cfg2 = config.clone();
        cfg2.band_width = 0;
        cfg2.adaptive_band = false;
        let mut graph2 = PoaGraph::new(reads[seed_idx], cfg2)?;
        for (i, read) in reads.iter().enumerate() {
            if i != seed_idx {
                graph2.add_read(read)?;
            }
        }
        return Ok(graph2);
    }

    Ok(graph)
}

/// Translate the `read_indices` on each [`Consensus`] from graph-internal
/// indices back to indices into the caller's input `reads` slice.
///
/// [`build_graph`] builds the graph as `PoaGraph::new(reads[seed_idx])`
/// (making the seed graph-internal index 0) followed by `add_read` on every
/// other read in input order. The internal-to-input mapping is therefore
/// `perm = [seed_idx] ++ (input indices 0..n, excluding seed_idx)`, i.e.
/// `input_idx = perm[internal_idx]`. This must mirror `build_graph`'s ordering
/// exactly; if that ordering ever changes, update this in lockstep.
///
/// Single-allele consensuses carry an empty `read_indices`, so remapping them
/// is a no-op. Only the free-function wrappers (which own the input slice and
/// the `seed_idx`) apply this; the stateful `PoaGraph::consensus_multi` leaves
/// `read_indices` in add-order, which is already the caller's own ordering.
fn remap_read_indices_to_input(consensuses: &mut [Consensus], seed_idx: usize, n_reads: usize) {
    let mut perm = Vec::with_capacity(n_reads);
    perm.push(seed_idx);
    for i in 0..n_reads {
        if i != seed_idx {
            perm.push(i);
        }
    }
    for c in consensuses.iter_mut() {
        for idx in c.read_indices.iter_mut() {
            *idx = perm[*idx];
        }
    }
}

/// Index of the read whose length is closest to the population's median
/// (by position in the length-sorted order, not interpolated) -- used by
/// `consensus_adaptive`'s seed-sensitivity retry to pick a re-seeding
/// candidate that is unlikely to be an atypically short (or long) outlier.
fn median_length_read_index(reads: &[&[u8]]) -> Option<usize> {
    if reads.is_empty() {
        return None;
    }
    let mut idx: Vec<usize> = (0..reads.len()).collect();
    idx.sort_by_key(|&i| reads[i].len());
    Some(idx[idx.len() / 2])
}

/// Shared implementation behind `consensus_adaptive`'s
/// `single_support_fraction > 0.3` branch and the standalone
/// [`consensus_fit_scored`]. Builds several candidate remedies from `reads`,
/// scores each against the actual read population with
/// [`analysis::consensus_fit`], and returns the best-fitting one. `c1` is the
/// caller's already-built pass-1 consensus (on `seed_idx`), passed in rather
/// than rebuilt so `consensus_adaptive` doesn't pay for a redundant graph
/// build it already has the result of.
///
/// See `consensus_adaptive`'s "Seed-sensitivity retry" doc section for the
/// full rationale and empirical validation notes.
fn seed_sensitivity_retry(
    reads: &[&[u8]],
    seed_idx: usize,
    config: &PoaConfig,
    c1: &Consensus,
) -> Result<AdaptiveResult, PoaError> {
    let mut candidates: Vec<(Consensus, AdaptiveAction)> =
        vec![(c1.clone(), AdaptiveAction::PassThrough)];

    // Remedy 1 (pre-existing): tighten the coverage floor, same seed.
    let mut cfg_tight = config.clone();
    cfg_tight.min_coverage_fraction = cfg_tight.min_coverage_fraction.max(0.6);
    if let Ok(c) = build_graph(reads, seed_idx, cfg_tight).and_then(|g| g.consensus()) {
        candidates.push((c, AdaptiveAction::NoisyTighten));
    }

    // Remedy 2: re-seed on a read near the population's median length
    // instead of whatever `seed_idx` the caller (or `select_seed::Auto`)
    // picked -- avoids the short-seed insert-fragmentation mechanism above
    // outright, rather than trying to recover from it after the fact.
    if let Some(median_idx) = median_length_read_index(reads) {
        if median_idx != seed_idx {
            if let Ok(c) =
                build_graph(reads, median_idx, config.clone()).and_then(|g| g.consensus())
            {
                candidates.push((c, AdaptiveAction::AlternateSeedRetry));
            }
        }
    }

    // Remedy 3: MajorityFrequency mode, same seed -- counts Delete
    // traversals explicitly instead of relying on edge weights, which
    // sometimes recovers length HeaviestPath loses to fragmentation (and
    // sometimes does not; that is exactly why this is scored rather than
    // assumed).
    let mut cfg_mf = config.clone();
    cfg_mf.consensus_mode = ConsensusMode::MajorityFrequency;
    if let Ok(c) = build_graph(reads, seed_idx, cfg_mf).and_then(|g| g.consensus()) {
        candidates.push((c, AdaptiveAction::MajorityFrequencyRetry));
    }

    // Score every candidate against the actual read population and keep the
    // best (lowest) fit score. Ties keep the earliest candidate, i.e. prefer
    // pass-1 over a rebuild, and prefer the pre-existing NoisyTighten remedy
    // over the two new ones, when scores tie exactly.
    let mut best = 0usize;
    let mut best_score = f64::INFINITY;
    for (i, (c, _)) in candidates.iter().enumerate() {
        let score = crate::analysis::consensus_fit(&c.sequence, reads, config);
        if score < best_score {
            best_score = score;
            best = i;
        }
    }
    let (chosen, action) = candidates.into_iter().nth(best).expect("non-empty");
    Ok(AdaptiveResult {
        consensuses: vec![chosen],
        action,
    })
}

/// Standalone seed-sensitivity check-and-retry, without `consensus_adaptive`'s
/// other passes (multi-allele bubble detection, truncation retry, semi-global
/// fallback).
///
/// Builds a pass-1 consensus exactly like [`consensus`]; if
/// `GraphStats::single_support_fraction` exceeds 0.3, additionally builds a
/// small set of candidate remedies and returns whichever one an empirical fit
/// score ([`analysis::consensus_fit`]) says the actual read population best
/// supports (see `consensus_adaptive`'s "Seed-sensitivity retry" doc section
/// for the full mechanism, rationale, and validated limitations). Otherwise
/// returns the pass-1 result unchanged with `action: PassThrough`.
///
/// Exists as a separate entry point (rather than requiring callers to go
/// through `consensus_adaptive`) for callers that want this specific retry
/// without `consensus_adaptive`'s other behavior -- e.g. this crate's own CLI,
/// which intentionally avoids `consensus_adaptive`'s multi-allele bubble
/// detection firing unexpectedly on a het SNP in a caller that expects
/// single-allele output.
pub fn consensus_fit_scored(
    reads: &[&[u8]],
    seed_idx: usize,
    config: &PoaConfig,
) -> Result<AdaptiveResult, PoaError> {
    validate(reads, seed_idx)?;
    let graph = build_graph(reads, seed_idx, config.clone())?;
    let stats = graph.stats();
    let c1 = graph.consensus()?;
    if stats.single_support_fraction > 0.3 {
        seed_sensitivity_retry(reads, seed_idx, config, &c1)
    } else {
        Ok(AdaptiveResult {
            consensuses: vec![c1],
            action: AdaptiveAction::PassThrough,
        })
    }
}

fn validate(reads: &[&[u8]], seed_idx: usize) -> Result<(), PoaError> {
    if reads.is_empty() {
        return Err(PoaError::EmptyInput);
    }
    if seed_idx >= reads.len() {
        return Err(PoaError::SeedOutOfBounds {
            index: seed_idx,
            len: reads.len(),
        });
    }
    Ok(())
}

// ── Public convenience wrappers ───────────────────────────────────────────────

/// Build a single-allele consensus from `reads`.
///
/// `seed_idx` is the index of the read used to initialise the graph; choose a
/// median-length read for best results.
pub fn consensus(
    reads: &[&[u8]],
    seed_idx: usize,
    config: &PoaConfig,
) -> Result<Consensus, PoaError> {
    validate(reads, seed_idx)?;
    build_graph(reads, seed_idx, config.clone())?.consensus()
}

/// Build a multi-allele consensus from `reads`.
///
/// Returns one [`Consensus`] per detected allele.  If no heterozygous bubble is
/// found the result is a single-element `Vec` equivalent to calling [`consensus`].
pub fn consensus_multi(
    reads: &[&[u8]],
    seed_idx: usize,
    config: &PoaConfig,
) -> Result<Vec<Consensus>, PoaError> {
    validate(reads, seed_idx)?;
    let mut alleles = build_graph(reads, seed_idx, config.clone())?.consensus_multi()?;
    remap_read_indices_to_input(&mut alleles, seed_idx, reads.len());
    Ok(alleles)
}

/// Two-pass adaptive consensus.
///
/// **Pass 1** builds a graph with the supplied config and computes [`GraphStats`].
/// **Pass 2** is selected by inspecting those stats:
///
/// | Condition | Pass-2 action |
/// |---|---|
/// | 1-3 bubbles, minority arm ≥ `min_allele_freq × n` | `consensus_multi` on pass-1 graph |
/// | Consensus < 60% of median input read length (banded, median ≤ 5000 bp) | Retry with `band_width = 0` |
/// | `single_support_fraction > 0.3` | Build several candidate remedies, keep whichever scores best against the actual reads (see below) |
/// | Coverage CV > 1.5 and mode is `Global` | Switch to `SemiGlobal`, rebuild |
/// | Otherwise | Return pass-1 single consensus; no rebuild |
///
/// Returns an [`AdaptiveResult`] containing `consensuses` (one element for
/// single-allele outcomes, two for diploid) and `action` (which pass-2 branch
/// fired, if any).  Inspect `action` to distinguish a clean pass-through from
/// a corrected result without re-running [`diagnose`].
///
/// ## Seed-sensitivity retry (`single_support_fraction > 0.3`)
///
/// This trigger fires on most genuinely repetitive/homogeneous graphs
/// (confirmed empirically to sit in the 0.3-0.45 range for the large
/// majority of periodic-repeat test scenarios, whether or not anything is
/// actually wrong), so it cannot be treated as "this specific consensus is
/// wrong" on its own -- it only means "this locus looks noisy/ambiguous
/// enough to be worth double-checking." Four candidates are built from the
/// *same* reads -- the untouched pass-1 result, a coverage-floor-tightened
/// rebuild (the sole pre-existing remedy), a rebuild re-seeded on a read
/// near the population's median length, and a
/// [`ConsensusMode::MajorityFrequency`] rebuild -- and scored against each
/// other with [`analysis::consensus_fit`] (mean per-read Insert+Delete ops
/// against each candidate). The lowest-scoring (best-fit) candidate is
/// returned, and `action` records which one won.
///
/// This is not a complete fix for the underlying seed-length sensitivity:
/// empirical validation across CAG/GAA repeats at several lengths, depths,
/// and error models found it reliably avoids ever making an
/// already-passing case *worse* (0 regressions across 21 tested passing
/// scenarios), and outright corrects some -- but not all -- confirmed
/// failure instances. On the residual failures, the chosen candidate is
/// sometimes no better than pass-1, and in at least one confirmed case
/// modestly *more* wrong than pass-1 by unit count (though that case was
/// already substantially wrong before this retry existed, not a regression
/// from a previously-correct result). No single absolute graph statistic
/// tested during development reliably distinguished the fixable cases from
/// the unfixable ones or from ordinary noise; see CHANGELOG for the full
/// empirical account, including the negative results.
///
/// ## Truncation detection
///
/// When banded DP is used on highly repetitive sequence (e.g. the AAAAG 5-mer
/// in RFC1), multiple DP diagonals score identically.  The band can lock onto a
/// shorter diagonal without approaching the band edge, producing a truncated
/// consensus with no error.  `consensus_adaptive` detects this by comparing the
/// pass-1 consensus length to the median input read length.  If the ratio is
/// below 0.60 and the median read length is at most 5000 bp, it retries with
/// unbanded alignment (`band_width = 0`), which forces the traceback to reach
/// the correct endpoint.  `band_width = 0` is genuinely O(read_len ×
/// graph_len) DP (see `PoaConfig::band_width`), so the retry is capped to
/// 5000 bp median read length -- the same cap the CLI's own truncation retry
/// uses -- rather than unconditionally retrying on arbitrarily long reads.
pub fn consensus_adaptive(
    reads: &[&[u8]],
    seed_idx: usize,
    config: &PoaConfig,
) -> Result<AdaptiveResult, PoaError> {
    validate(reads, seed_idx)?;

    // ── Pass 1 ───────────────────────────────────────────────────────────────
    let graph = build_graph(reads, seed_idx, config.clone())?;
    let stats = graph.stats();

    // ── Decision ─────────────────────────────────────────────────────────────
    let n = reads.len();
    let allele_threshold = (n as f64 * config.min_allele_freq).ceil() as usize;

    // Multi-allele: few bubbles with a well-supported minority arm.
    // Re-use the pass-1 graph; consensus_multi rebuilds per-allele sub-graphs.
    if stats.bubble_count >= 1
        && stats.bubble_count <= 3
        && stats.max_bubble_depth >= allele_threshold
    {
        let mut consensuses = graph.consensus_multi()?;
        remap_read_indices_to_input(&mut consensuses, seed_idx, n);
        return Ok(AdaptiveResult {
            consensuses,
            action: AdaptiveAction::MultiAllele,
        });
    }

    // Compute pass-1 consensus; used for truncation check and as the fallthrough
    // return value if no second pass fires.
    let c1 = graph.consensus()?;

    // Truncation: consensus much shorter than median read length.
    // Indicates banded DP converged to the wrong diagonal (bug #4 pattern:
    // degenerate scoring landscape in highly repetitive sequence means the band
    // can snap to an off-by-N×period diagonal without approaching the edge).
    // Unbanded forces the traceback to terminate at the correct graph endpoint.
    //
    // band_width=0 (paired with adaptive_band=false below) is genuinely
    // unbanded O(read_len * graph_len) DP, so the retry is gated on read length
    // the same way the CLI gates its own truncation retry: above 5000 bp the
    // memory/time cost is no longer a safe unconditional fallback, so this
    // returns the (still truncated) pass-1 result rather than retrying.
    let was_banded = config.band_width > 0 || config.adaptive_band;
    if was_banded {
        let median_len = c1.graph_stats.median_input_read_len;
        if median_len > 0
            && median_len <= 5_000
            && (c1.sequence.len() as f64) < 0.6 * median_len as f64
        {
            let mut cfg2 = config.clone();
            cfg2.band_width = 0;
            cfg2.adaptive_band = false;
            let c2 = build_graph(reads, seed_idx, cfg2)?.consensus()?;
            let recovered = (c2.sequence.len() as f64) >= 0.6 * median_len as f64;
            return Ok(AdaptiveResult {
                consensuses: vec![c2],
                action: AdaptiveAction::TruncationRetry { recovered },
            });
        }
    }

    // Noisy input: high fraction of singleton-supported nodes.
    //
    // Before this pass existed, the sole remedy was tightening the coverage
    // floor and rebuilding on the *same* seed. That does not address a
    // different, confirmed root cause behind this same trigger: an
    // auto-selected seed that happens to be atypically short relative to the
    // true read population can make heaviest_path/the interior filter
    // systematically under-call a periodic/homogeneous repeat's length --
    // the extra content most reads carry over the short seed has to be
    // inserted somewhere, and in a homogeneous repeat any position is an
    // equally valid place to insert it, so different reads scatter their
    // inserts and no single insertion position accumulates enough coverage
    // to survive on its own (confirmed via ground-truth read-length
    // comparison on synthetic CAG/GAA repeats; reproduces identically fully
    // unbanded, so it is not a band-width artifact). Tightening the coverage
    // floor cannot fix this: it filters harder against evidence that is
    // already fragmented, it does not un-fragment it.
    //
    // No single absolute, graph-level statistic reliably distinguished a
    // genuinely-too-short consensus from ordinary sequencing-noise-heavy
    // ones across scenario families in empirical testing (tried:
    // `single_support_fraction`, `bubble_count`, `edge_weight_gini`,
    // `mean_column_entropy`, seed length relative to the read population's
    // median/longest read, and an aggregate "off-spine fork mass" measure --
    // none separated cleanly; see CHANGELOG for the full account). What DID
    // discriminate, empirically, is comparing candidate consensuses against
    // each other using `consensus_fit` (mean per-read Insert+Delete ops
    // against a candidate, normalised by length) -- i.e. scoring, not
    // guessing, and always keeping whichever candidate (including the
    // untouched pass-1 result) the actual read population best supports.
    if stats.single_support_fraction > 0.3 {
        return seed_sensitivity_retry(reads, seed_idx, config, &c1);
    }

    // Uneven boundary coverage (high CV) suggests partial reads in global mode.
    // Switch to semi-global alignment and rebuild.
    let cv = if stats.coverage_mean > 0.0 {
        stats.coverage_variance.sqrt() / stats.coverage_mean
    } else {
        0.0
    };
    if cv > 1.5 && config.alignment_mode == AlignmentMode::Global {
        let mut cfg2 = config.clone();
        cfg2.alignment_mode = AlignmentMode::SemiGlobal;
        return Ok(AdaptiveResult {
            consensuses: vec![build_graph(reads, seed_idx, cfg2)?.consensus()?],
            action: AdaptiveAction::SemiGlobalFallback,
        });
    }

    // No second pass needed.
    Ok(AdaptiveResult {
        consensuses: vec![c1],
        action: AdaptiveAction::PassThrough,
    })
}

/// Build a consensus from two non-overlapping read groups with a gap of
/// unknown length between them.
///
/// Use this when reads cover only the left end and right end of a template
/// with no read bridging the middle — for example, when a repeat expansion
/// is longer than any individual read.
///
/// Each group is assembled independently using [`consensus`].  The results
/// are concatenated and a single [`CoverageGap`] with
/// [`GapKind::Unknown`] is inserted at the join point.
///
/// The returned `Consensus.n_reads` is the sum of both groups.
/// `Consensus.graph_stats` reflects the left group's graph.
///
/// # Size interpretation
///
/// ```text
/// |=left_consensus=|???unknown???|=right_consensus=|
/// total ≥ left_consensus.len() + right_consensus.len()
/// ```
///
/// `gap.min_size()` returns `None` (size unknown).  The caller can report:
/// ```text
/// format!("≥{} bp with a gap of unknown length",
///     cons.sequence.len())
/// ```
pub fn bridged_consensus(
    left_reads: &[&[u8]],
    left_seed_idx: usize,
    right_reads: &[&[u8]],
    right_seed_idx: usize,
    config: &PoaConfig,
) -> Result<Consensus, PoaError> {
    validate(left_reads, left_seed_idx)?;
    validate(right_reads, right_seed_idx)?;

    let left = build_graph(left_reads, left_seed_idx, config.clone())?.consensus()?;
    let right = build_graph(right_reads, right_seed_idx, config.clone())?.consensus()?;

    let join = left.sequence.len();

    let mut sequence = left.sequence.clone();
    sequence.extend_from_slice(&right.sequence);

    let mut coverage = left.coverage.clone();
    coverage.extend_from_slice(&right.coverage);

    let mut path_weights = left.path_weights.clone();
    path_weights.extend_from_slice(&right.path_weights);

    // Carry forward any spanning gaps from each segment (offset right gaps by join).
    let mut gaps = left.gaps.clone();
    for gap in &right.gaps {
        gaps.push(CoverageGap {
            start: gap.start + join,
            end: gap.end + join,
            kind: gap.kind,
        });
    }
    // Insert the unknown gap at the join point, then sort by position.
    gaps.push(CoverageGap {
        start: join,
        end: join,
        kind: GapKind::Unknown,
    });
    gaps.sort_by_key(|g| g.start);

    Ok(Consensus {
        sequence,
        coverage,
        path_weights,
        n_reads: left.n_reads + right.n_reads,
        graph_stats: left.graph_stats,
        gaps,
        bubble_sites: vec![],
        read_indices: vec![],
    })
}