use crate::graph::AlignOp;
use crate::{self as poa_consensus, AlignmentMode, ConsensusMode, PoaConfig, PoaError, PoaGraph};
fn b(s: &str) -> Vec<u8> {
s.as_bytes().to_vec()
}
fn s(v: &[u8]) -> String {
String::from_utf8_lossy(v).into_owned()
}
/// Build consensus from reads; seed_idx selects the first read added to the graph.
fn consensus(reads: &[Vec<u8>], seed_idx: usize) -> Vec<u8> {
let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
for (i, read) in reads.iter().enumerate() {
if i == seed_idx {
continue;
}
graph.add_read(read).unwrap();
}
graph.consensus().unwrap().sequence
}
fn consensus_cfg(reads: &[Vec<u8>], seed_idx: usize, cfg: PoaConfig) -> Vec<u8> {
let mut graph = PoaGraph::new(&reads[seed_idx], cfg).unwrap();
for (i, read) in reads.iter().enumerate() {
if i == seed_idx {
continue;
}
graph.add_read(read).unwrap();
}
graph.consensus().unwrap().sequence
}
// ── Error cases ───────────────────────────────────────────────────────────────
#[test]
fn empty_reads() {
let result = PoaGraph::new(&[], PoaConfig::default());
assert!(
matches!(result, Err(PoaError::EmptyInput)),
"expected EmptyInput"
);
}
#[test]
fn below_min_reads() {
let cfg = PoaConfig {
min_reads: 3,
..Default::default()
};
let mut graph = PoaGraph::new(&b("ACGT"), cfg).unwrap();
graph.add_read(&b("ACGT")).unwrap();
let result = graph.consensus();
assert!(
matches!(result, Err(PoaError::InsufficientDepth { got: 2, min: 3 })),
"expected InsufficientDepth, got {:?}",
result
);
}
#[test]
fn seed_out_of_bounds() {
// Callers must check seed_idx < reads.len() before calling PoaGraph::new.
// Document the expected guard pattern.
let reads = vec![b("ACGT"), b("ACGT")];
let idx = 5usize;
assert!(idx >= reads.len(), "caller must guard seed_idx before use");
}
// ── Basic correctness ─────────────────────────────────────────────────────────
#[test]
fn single_read_passthrough() {
let reads = vec![b("CATCATCAT")];
assert_eq!(consensus(&reads, 0), b("CATCATCAT"));
}
#[test]
fn two_identical_reads() {
let reads = vec![b("CATCATCAT"), b("CATCATCAT")];
assert_eq!(consensus(&reads, 0), b("CATCATCAT"));
}
#[test]
fn majority_base_wins() {
let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CGTCATCAT")];
assert_eq!(s(&consensus(&reads, 0)), "CATCATCAT");
}
#[test]
fn single_outlier_not_inflated() {
let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCATCAT")];
assert_eq!(consensus(&reads, 0).len(), 9);
}
#[test]
fn length_variation_longer_wins() {
let reads = vec![b("CATCATCATCAT"), b("CATCATCATCAT"), b("CATCATCAT")];
assert_eq!(consensus(&reads, 0).len(), 12);
}
#[test]
fn no_inflation_with_length_noise() {
let reads = vec![
b("CAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAG"),
];
assert_eq!(consensus(&reads, 0).len(), 15);
}
#[test]
fn no_inflation_phox2b_like() {
let reads = vec![
b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
b("GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA"),
];
assert_eq!(consensus(&reads, 0).len(), 60);
}
#[test]
fn single_base_reads() {
let reads = vec![b("A"), b("A"), b("A")];
assert_eq!(consensus(&reads, 0), b("A"));
}
// ── Boundary trim ─────────────────────────────────────────────────────────────
#[test]
fn boundary_trim_leading_seed_artifact() {
let reads = vec![
b("XXXCATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "CATCATCAT", "got: {}", result);
}
#[test]
fn boundary_trim_trailing_seed_artifact() {
let reads = vec![
b("CATCATCATXXX"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "CATCATCAT", "got: {}", result);
}
// ── Diagnostic tests from ref/poa.rs ─────────────────────────────────────────
#[test]
fn diag_sca3_t3_tail_seed_t1() {
let reads = vec![
b("CAGCAGCAGT"),
b("CAGCAGCAGTTT"),
b("CAGCAGCAGTTT"),
b("CAGCAGCAGTTT"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "CAGCAGCAGTTT", "got: {}", result);
}
#[test]
fn diag_sca3_t3_tail_seed_t3() {
let reads = vec![
b("CAGCAGCAGTTT"),
b("CAGCAGCAGTTT"),
b("CAGCAGCAGTTT"),
b("CAGCAGCAGT"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "CAGCAGCAGTTT", "got: {}", result);
}
#[test]
fn diag_sca31_trailing_interrupt_seed_missing() {
let reads = vec![
b("ATTATTATTATT"),
b("ATTATTATTATTATA"),
b("ATTATTATTATTATA"),
b("ATTATTATTATTATA"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "ATTATTATTATTATA", "got: {}", result);
}
#[test]
fn diag_sca3_interrupt_position_single_outlier() {
let maj = b("CAGCAGCAGCAGCAGGTTCAGCAG");
let out = b("CAGCAGCAGCAGCAGCAGGTTCAGCAG");
let reads = vec![maj.clone(), maj.clone(), maj.clone(), out];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "CAGCAGCAGCAGCAGGTTCAGCAG", "got: {}", result);
}
#[test]
fn diag_sca8_minority_trailing_extension_trimmed() {
let base = b("CAGCAGCAGCAGCAG");
let extend = b("CAGCAGCAGCAGCAGGCT");
let reads = vec![
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
extend.clone(),
extend.clone(),
extend.clone(),
];
assert_eq!(consensus(&reads, 0).len(), 15);
}
#[test]
fn diag_sca8_minority_trailing_extension_seed_extends() {
let base = b("CAGCAGCAGCAGCAG");
let extend = b("CAGCAGCAGCAGCAGGCT");
let reads = vec![
extend.clone(),
extend.clone(),
extend.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
base.clone(),
];
assert_eq!(consensus(&reads, 0).len(), 15);
}
#[test]
fn diag_sca31_trailing_interrupt_before_flank() {
let flank = b("GCGCGCGC");
let mut seed_read = b("ATTATTATTATT");
seed_read.extend_from_slice(&flank);
let mut maj_read = b("ATTATTATTATTATA");
maj_read.extend_from_slice(&flank);
let reads = vec![
seed_read,
maj_read.clone(),
maj_read.clone(),
maj_read.clone(),
];
let result = s(&consensus(&reads, 0));
let expected: String = "ATTATTATTATTATA"
.chars()
.chain("GCGCGCGC".chars())
.collect();
assert_eq!(result, expected, "got: {}", result);
}
#[test]
fn diag_sca3_interrupt_position_long_repeat_with_flank() {
let flank = b("CTGCTGCTG");
let make = |repeat_pre: &str, interrupt: &str, repeat_post: &str| -> Vec<u8> {
let mut v = repeat_pre.as_bytes().to_vec();
v.extend_from_slice(interrupt.as_bytes());
v.extend_from_slice(repeat_post.as_bytes());
v.extend_from_slice(&flank);
v
};
let maj = make("CAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAGCAG");
let out = make("CAGCAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAG");
let reads = vec![maj.clone(), maj.clone(), maj.clone(), out];
let result = s(&consensus(&reads, 0));
let expected = s(&make("CAGCAGCAGCAGCAGCAGCAGCAG", "GTT", "CAGCAGCAG"));
assert_eq!(result, expected, "got: {}", result);
}
#[test]
#[ignore]
fn diag_frda_gaa_rotation_phase() {
let gaa_phase = b("GAAGAAGAAGAA");
let aag_phase = b("AAGAAGAAGAAG");
let aga_phase = b("AGAAGAAGAAGA");
let reads = vec![gaa_phase.clone(), gaa_phase.clone(), aag_phase, aga_phase];
let result = s(&consensus(&reads, 0));
assert_eq!(result.len(), 12, "got: '{}'", result);
}
#[test]
fn diag_frda_gaa_rotation_with_flanking() {
let make = |repeat: &str| -> Vec<u8> {
let mut v = b("TTTCCC");
v.extend_from_slice(repeat.as_bytes());
v.extend_from_slice(b("GGGAAA").as_slice());
v
};
let reads = vec![
make("GAAGAAGAAGAA"),
make("GAAGAAGAAGAA"),
make("AAGAAGAAGAAG"),
make("AGAAGAAGAAGA"),
];
// All reads are constructed to span the full TTTCCC..GGGAAA region, so global
// alignment is correct here: the traceback is forced to (n, l), anchoring both
// ends of the consensus. With SemiGlobal the traceback exits from the max score
// in the last row; rotation-phase bubbles create near-tied scores at the tail and
// the traceback exits one base early. Real data uses SemiGlobal because reads do
// not always span the full locus; here they do, so global is unambiguous.
let cfg = PoaConfig {
band_width: 0,
adaptive_band: false,
alignment_mode: AlignmentMode::Global,
..PoaConfig::default()
};
let result = s(&consensus_cfg(&reads, 0, cfg));
assert_eq!(result.len(), 24, "got: '{}'", result);
}
#[test]
fn diag_phase_shift_first_node_coverage() {
let reads = vec![
b("GAAGAA"),
b("GAAGAA"),
b("GAAGAA"),
b("GAAGAA"),
b("AAGAAG"),
];
let result = s(&consensus(&reads, 0));
assert_eq!(result, "GAAGAA", "got: '{}'", result);
}
#[test]
fn diag_phase_shift_majority_trims_first_base() {
let reads = vec![
b("GAAGAA"),
b("GAAGAA"),
b("AAGAAG"),
b("AAGAAG"),
b("AAGAAG"),
];
let result = s(&consensus(&reads, 0));
// The majority is the same sequence in a different phase; correct length is 6.
assert_eq!(result.len(), 6, "got: '{}'", result);
}
#[test]
fn diag_dab1_sca37_attttc_lookahead_arm_length_bias() {
// DAB1 SCA37 hap2, HiFi sample, chr1:57367043-57367121 +/- 20bp flank
// (28 reads). Read index 0 has a genuinely different repeat-unit count
// (8x AAAAT vs the majority's 7x), which creates a 2-arm bubble at the
// seed backbone: a 1-node "no insertion" arm vs read 0's ~19-node
// insertion detour. The lookahead resolver used to gate on whether the
// *longest* arm could provide LOOKAHEAD_K bases, then score each arm only
// over its own available length -- so the 1-node arm could score at most
// 1 x match_score while the 19-node arm scored over the full
// LOOKAHEAD_K bases, guaranteeing the long arm always won regardless of
// which one actually matched. Every subsequent read landed on the same
// over-scored arm, fragmenting read support and silently dropping bases
// scattered through the AAAAT run. The bug was highly seed-dependent:
// some seeds happened to avoid read 0's bubble contamination and produced
// the correct consensus, others did not. This test pins seed=1 (which
// produced ACAGTGAGACCCTGTCTCCAAAAATAAAATAAAATAAAATAAAAAAAAAAATAAATTAGCCGCATGGT
// before the fix, a 68bp corrupted consensus dropping 10bp across three
// sites) as a regression guard.
let reads = vec![
b("ACAGTGAGACCCTGTCTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATTAGCCCAGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAATAAATAAATAAATAAATAAATAAATAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ATCGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b(
"ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT",
),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATTAGCCAGGCATGGT"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAATAAATTA"),
b("ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT"),
];
let expected = s(&b(
"ACAGTGAGACCCTGTCTCCACAAAATAAAATAAAATAAAATAAAATAAAATAAAATAAATAAATTAGCCAGGCATGGT",
));
for seed_idx in 0..reads.len() {
let result = s(&consensus(&reads, seed_idx));
assert_eq!(result, expected, "seed={} got: {}", seed_idx, result);
}
}
#[test]
fn diag_dmd_ctt_repeat_interior_filter_global_threshold_bias() {
// DMD hap2, HG02968 HiFi, chrX:31284557-31284613 +/- 20bp flank
// (15 reads, extracted via bedpull). This is the CTT-repeat counterpart
// to the DAB1 test above, but the root cause is in the interior filter,
// not the aligner. At seed=12, a node on heaviest_path's own winning
// route has real local majority support at its own bubble, but its
// Match coverage falls short of the *global* min_cov purely because
// other reads diverged at an earlier bubble. Dropping it spliced its
// flanking nodes together, fabricating a junction no read ever
// produced (an "AAACTGCAATAACGA" 5' boundary that no read's actual
// flank matches). See `diag_dmd_ctt_majority_delete_residual` below for
// a related, still-open case this fix does not cover.
let reads = vec![
b(
"AATAATTAAATACTGTTTTTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTTTTTTGTAGAGGTGGGGTCT",
),
b(
"AGAACTGCAATAACGAACTGTCTCTCTTTCTTCTTCTTCTCCTCCCTCCCTCCCTCCTCTCTTCTCTCTTCCTCCCTCCTCCCTCCCTCCCTCCCTCCCTCCTCTCCTCCCTCCTCTCCTCTCTTCTCTTCTTTCTTCTTCTCTTCTCTTCTCCTCTCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCTTTTTTTTTTTTTGGCAGAGGTGGTGTCT",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"ACAACTGCAATAACCAGACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAGAACGAACTGTTTCCTCCTTCTCCTCCTCCTTCTTCTTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTTCTTCTTTGGCAGAGGTGGGGGTGTTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCCTCTTCTTCTTCTTCTTCTTCTTCTTCTCTTCTTCTTCTCTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTTAC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTTTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTCAGAGGTGACATGT",
),
b(
"AGACTGCAATAAGGGACTTCTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTCTTTCTTTTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTCTTCTTCTTTCTTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTTTTTGGCGAGGTGGAGTGC",
),
b(
"AGAAACTGCAGAACGACTGTTTCTTCTTCTTCTTCTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTCCCTCTCTTTTTGGCAGAGGTGGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
];
let cfg = PoaConfig {
min_coverage_fraction: 0.6,
..PoaConfig::default()
};
let result = s(&consensus_cfg(&reads, 12, cfg));
assert!(
!result.contains("CTTTTC") && !result.starts_with("AAACTGCAATAACGA"),
"interior filter fabricated sequence not present in any read: {}",
result
);
}
#[test]
// Formerly #[ignore]d as a known residual of the majority-Delete interior-filter
// fabrication (Known Bug #6). Now a live regression test: the fabrication is
// genuinely fixed. Bisected — this exact read set + seed produced the impossible
// "CTTTTC" run at commit d24f133 (test's own introduction) and still at 19b81e8,
// then went green at 64a839f ("RFC1 leading interrupt bug", 2026-07-08), the
// backward-fork-search rescue (Known Bug #8, FORK_SEARCH_HOPS=64). Fixed by that
// earlier interior-filter work, NOT by the later bypass-edge Delete rework.
fn diag_dmd_ctt_majority_delete_residual() {
let reads = vec![
b(
"AATAATTAAATACTGTTTTTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTTTTTTGTAGAGGTGGGGTCT",
),
b(
"AGAACTGCAATAACGAACTGTCTCTCTTTCTTCTTCTTCTCCTCCCTCCCTCCCTCCTCTCTTCTCTCTTCCTCCCTCCTCCCTCCCTCCCTCCCTCCCTCCTCTCCTCCCTCCTCTCCTCTCTTCTCTTCTTTCTTCTTCTCTTCTCTTCTCCTCTCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCCCTCTTTTTTTTTTTTTGGCAGAGGTGGTGTCT",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"ACAACTGCAATAACCAGACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAGAACGAACTGTTTCCTCCTTCTCCTCCTCCTTCTTCTTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTCCTTCTTCTTTGGCAGAGGTGGGGGTGTTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCCTCTTCTTCTTCTTCTTCTTCTTCTTCTCTTCTTCTTCTCTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTTAC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTTTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTCAGAGGTGACATGT",
),
b(
"AGACTGCAATAAGGGACTTCTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTCTTTCTTTTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTCTTCTTCTTTCTTTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTCTTTTTTTGGCGAGGTGGAGTGC",
),
b(
"AGAAACTGCAGAACGACTGTTTCTTCTTCTTCTTCTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTCCCTCTCTTTTTGGCAGAGGTGGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
b(
"AGAACTGCAATAACGAACTGTTTTTTCTTCTTCTTCTTCTTTCTTCTTCCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTCTTTTTTTTGGCAGAGGTGGGGTGTCTC",
),
];
let cfg = PoaConfig {
min_coverage_fraction: 0.6,
..PoaConfig::default()
};
let result = s(&consensus_cfg(&reads, 1, cfg));
assert!(
!result.contains("CTTTTC"),
"interior filter fabricated an impossible CTT-repeat run: {}",
result
);
}
#[test]
fn diag_rfc1_aaaag_spurious_g_interrupt_local_rescue_noise() {
// CANVAS_RFC1 (AAAAG pentanucleotide repeat), Hap1, HG002 real data (10
// padded, medoid-filtered reads at seed=6, as bladerunner fed them to
// consensus_adaptive). Reported bug: consensus produced
// (AAAAG)30(G)1(AAAAG)86 -- a single extra "G" interrupt after unit 30
// that no read supports at that position (checked all 10 reads at the
// equivalent offset; only 2 of 10 have an "AAAAGG" anywhere at all, and
// both are far from unit 30).
//
// Root cause: the local-dominance rescue added for the DAB1/DMD
// interior-filter bug (see diag_dab1_sca37_attttc_lookahead_arm_length_bias
// and diag_dmd_ctt_repeat_interior_filter_global_threshold_bias above)
// gates only on relative local majority (a node's Match coverage clears
// a majority of its predecessor's *local* out-edge weight total), not on
// the absolute size of that local population. In this AAAAG region,
// repeated small phase-registration bubbles fragment the read pool
// rapidly; by the time this specific fork is reached, only 3 of the 10
// reads are even still on this exact sub-path. A 2-vs-1 split among
// those 3 clears the *local* majority bar trivially (local_min_cov=2)
// even though it is statistical noise from a fragmented region, not a
// real minority allele the global threshold unfairly suppresses. Fixed
// by additionally requiring the bubble's own local population to clear
// the *global* min_cov before trusting a local majority within it.
let reads = vec![
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAGAAAAGAAAAGAAGAAAGAAAGAAGAAAAGAAAGAAGAAAGAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAGGGAAAGAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAGAAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAAGAAAAGAAAGAAAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAGAAAAGAAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAA",
),
];
for seed_idx in 0..reads.len() {
let result = s(&consensus(&reads, seed_idx));
assert!(
!result.contains("AAAAGG"),
"seed={} produced an unsupported G interrupt: {}",
seed_idx,
result
);
}
}
#[test]
fn diag_rfc1_leading_interrupt_partial_read_population_accounting() {
// CANVAS_RFC1 (AAAAG pentanucleotide repeat), Hap2, chm13 real data
// (6 reads: 4 full-length spanning reads with independent single-base
// noise at scattered positions, plus 2 short partial reads covering
// only the first ~55-60bp). The only signal all 4 full reads agree on
// is a single "AAG" interrupt 12 clean units before the 3' boundary;
// everything else is per-read noise at different positions in each.
//
// Root cause: the interior filter's local-dominance rescue (added for
// the DAB1/DMD/RFC1-G-interrupt bugs above) only checked the
// *immediate* predecessor for a fork. Once the 2 short partial reads
// dropped out (their alignment simply ends there under semi-global
// mode -- a legitimate reduction in local population, not a bubble),
// `n_reads` stayed pinned at 6 for the *global* min_cov, even though
// only 4 reads structurally reach deeper positions. A node one hop
// downstream of an already-rescued fork (whose own predecessor has
// only one out-edge, i.e. isn't itself a fork) always fell through to
// the global check and got dropped, even when it was the clear
// majority (3 of the 4 still-active reads). Confirmed the 4-full-reads
// in isolation produce a perfectly clean consensus with no fabrication
// at all in this region.
//
// Fixed by walking backward (bounded) to find the nearest real fork
// rather than requiring the immediate predecessor to be one, so a
// rescued node's own single successor can inherit its fork's
// population instead of being judged against the global read count.
// A second fix in the same area (see
// diag_rfc1_leading_interrupt_delete_driven_forkless_gap below) covers
// the case where no fork exists at all anywhere nearby.
//
// This test covers only the *first ~230bp*, which the 4-full-reads
// diagnostic confirmed is fully explained by clean population
// accounting. Positions beyond that hit a different, deeper issue --
// genuine AAAAG periodic-alignment ambiguity (already documented as
// Known Bug #3/#4 in CLAUDE.md) -- tracked separately, not by this
// test.
let reads = vec![
b(
"AAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b("AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA"),
b("AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA"),
b(
"AAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGA",
),
];
let result = s(&consensus(&reads, 0));
let region = &result[..230.min(result.len())];
let bytes = region.as_bytes();
let g_positions: Vec<usize> = bytes
.iter()
.enumerate()
.filter(|&(_, &b)| b == b'G')
.map(|(i, _)| i)
.collect();
for w in g_positions.windows(2) {
let gap = w[1] - w[0] - 1;
assert_eq!(
gap, 4,
"unsupported gap of {} (expected 4) between G at {} and {} in: {}",
gap, w[0], w[1], region
);
}
}
#[test]
fn diag_rfc1_leading_interrupt_delete_driven_forkless_gap() {
// Same RFC1 AAAAG hap2 read set as the test above, but covering the
// *whole* consensus rather than just the first ~230bp.
//
// Root cause: one read (read index 2 here) individually Deletes a
// single "G" out of an otherwise fully unbranched, unanimous run of
// matches -- no fork is ever created, because Match and Delete share
// the same edge; only the node's own delete_count records that one
// read skipped it. The backward-fork-search rescue above only fires
// when a fork exists somewhere nearby; here there is none at all for
// dozens of nodes in either direction, so the node fell through to the
// *global* min_cov (based on all 6 reads, including the 2 short reads
// that never reach this deep) and was dropped even though 3 of the 4
// reads that actually reach this point matched it. Removing that one
// "G" merged the two flanking 4-A groups into a single 8-A run in the
// output -- reported as an "(A)4(AAAAG)1(AAAG)1" interrupt with no
// read support at that position.
//
// Fixed by falling back, when no fork is found, to the node's own
// coverage + delete_count as the local population (every read that
// reaches an unforked node either matches or deletes it, so this sum
// is exact) -- gated on that sum itself clearing the *global* min_cov,
// the same way the fork branch gates on the fork's total. That gate is
// what keeps this fix from also pulling in a genuine trailing/leading
// extension (a minority of reads simply longer than the rest, with
// nothing downstream to reconverge with -- see
// edge_extreme_length_variation_majority_wins in tests/sv_analysis.rs
// and long_repeat_length_majority_wins above, both regressed by an
// earlier version of this fix that used a forward-recovery scan
// instead, which a stray noisy read overlapping the tail could fool).
let reads = vec![
b(
"AAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAGAAAAGAAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b(
"AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA",
),
b("AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA"),
b("AAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAA"),
b(
"AAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGAAAAGA",
),
];
let result = s(&consensus(&reads, 0));
let bytes = result.as_bytes();
let g_positions: Vec<usize> = bytes
.iter()
.enumerate()
.filter(|&(_, &b)| b == b'G')
.map(|(i, _)| i)
.collect();
// The fixed fabrication produced gaps of 8 and 3 (from dropping one G
// out of two adjacent 4-A groups). Assert those specific patterns are
// gone. A single gap of 2 elsewhere is the known, separately-tracked
// periodic-alignment-ambiguity residual (see CLAUDE.md Known Bugs
// #3/#4) and is deliberately not asserted away here.
for w in g_positions.windows(2) {
let gap = w[1] - w[0] - 1;
assert!(
gap != 8 && gap != 3,
"delete-driven forkless gap fabrication reappeared: gap {} between G at {} and {} in: {}",
gap,
w[0],
w[1],
result
);
}
}
#[test]
fn diag_sca3_t3_tail_with_flank() {
let flank = b("CCTCCTCCT");
let make = |tail: &str| -> Vec<u8> {
let mut v = b("CAGCAGCAG");
v.extend_from_slice(tail.as_bytes());
v.extend_from_slice(&flank);
v
};
let reads = vec![make("T"), make("TTT"), make("TTT"), make("TTT")];
let result = s(&consensus(&reads, 0));
let expected = s(&make("TTT"));
assert_eq!(result, expected, "got: {}", result);
}
// ── Real-data reproduction ────────────────────────────────────────────────────
#[test]
fn diag_sca8_real_sequences_no_flank() {
let maj57 = b("TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT");
let min75 = b("TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT");
let min81 =
b("TACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT");
let mut reads: Vec<Vec<u8>> = std::iter::repeat(maj57.clone()).take(32).collect();
reads.push(b(
"TACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTAC",
));
reads.push(b(
"TACTACTACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCT",
));
reads.push(b("TACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT"));
reads.push(min75.clone());
reads.push(min81.clone());
reads.push(min81.clone());
let result = consensus(&reads, 0);
assert_eq!(
result.len(),
maj57.len(),
"SCA8 consensus must match majority length {}, got len {}: '{}'",
maj57.len(),
result.len(),
s(&result)
);
}
#[test]
fn diag_sca8_real_sequences_with_flank() {
let flank_l = b("GCTTCGAAGTC");
let flank_r = b("AAACGGTTCCA");
let make = |repeat: &[u8]| -> Vec<u8> {
let mut v = flank_l.clone();
v.extend_from_slice(repeat);
v.extend_from_slice(&flank_r);
v
};
let maj = make(&b(
"TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCT",
));
let min75 = make(&b(
"TACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT",
));
let min81 = make(&b(
"TACTACTACTACTACTACTACTACTACTACTACTACTACTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCT",
));
let mut reads: Vec<Vec<u8>> = std::iter::repeat(maj.clone()).take(32).collect();
reads.push(min75);
reads.push(min81.clone());
reads.push(min81);
let result_len = consensus(&reads, 0).len();
assert_eq!(
result_len,
maj.len(),
"SCA8 flanked: got {}, expected {}",
result_len,
maj.len()
);
}
// ── Banded DP ─────────────────────────────────────────────────────────────────
#[test]
fn banded_matches_unbanded_small() {
// Banded and unbanded must produce identical results when the band is wide
// enough to cover the optimal path.
let reads = vec![
b("CAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAGCAG"),
b("CAGCAGCAGCAGCAG"),
];
let unbanded = consensus(&reads, 0);
let cfg_banded = PoaConfig {
band_width: 50,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg_banded);
assert_eq!(unbanded, banded, "banded vs unbanded mismatch");
}
#[test]
fn adaptive_band_matches_unbanded() {
let reads = vec![
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCATCAT"),
b("CATCATCAT"),
];
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
adaptive_band: true,
adaptive_band_b: 5,
adaptive_band_f: 0.1,
..Default::default()
};
let adaptive = consensus_cfg(&reads, 0, cfg);
assert_eq!(unbanded, adaptive, "adaptive band vs unbanded mismatch");
}
#[test]
fn band_too_narrow_fallback_to_unbanded() {
// seed = 1 A, read = 30 A's: the 2-pass banded retry exhausts all banded
// options but the 3-pass unbanded fallback recovers. BandTooNarrow is now
// an internal signal, not a user-visible error.
let seed = b("A");
let read = b("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"); // 30 A's
let cfg = PoaConfig {
band_width: 2,
..Default::default()
};
let mut graph = PoaGraph::new(&seed, cfg).unwrap();
let result = graph.add_read(&read);
assert!(
result.is_ok(),
"3-pass retry must recover via unbanded fallback, got {:?}",
result.map(|_| ())
);
}
#[test]
fn large_length_variance_banded() {
// 3 reads of 15 bp majority + 1 read of 60 bp outlier.
// Band must be wide enough to cover the 45-base insertion from the outlier.
// With band_width=50 the banded result should match unbanded.
let maj = b("CAGCAGCAGCAGCAG");
let outlier = b("CAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAG");
let reads = vec![maj.clone(), maj.clone(), maj.clone(), outlier];
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 50,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded.len(),
unbanded.len(),
"banded length mismatch with large variance"
);
assert_eq!(
banded, unbanded,
"banded result mismatch with large variance"
);
}
// ── Semi-global alignment ─────────────────────────────────────────────────────
#[test]
fn partial_reads_semi_global() {
// 3 full reads + 1 partial: full reads are the majority so trailing region
// has coverage 3 >= min_cov=3, consensus is the full sequence.
let full = b("ACGTACGTACGT");
let partial = b("ACGTACGT");
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default()
};
let reads = vec![full.clone(), full.clone(), full.clone(), partial];
let result = consensus_cfg(&reads, 0, cfg);
assert_eq!(
result.len(),
12,
"semi-global: got len {}, seq: '{}'",
result.len(),
s(&result)
);
}
// ── Partial read coverage behaviour ──────────────────────────────────────────
#[test]
fn one_spanning_many_partial_default_min_cov_truncates() {
// With default min_cov (≈ n/2 + 1), a single spanning read never provides
// enough coverage to keep boundary nodes when partial reads dominate.
// This test documents the known behaviour so a future change doesn't
// silently alter it.
let spanning = b("ACGTACGTACGT"); // 12 bp
let partial = b("ACGTACGT"); // 8 bp (covers prefix only)
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default()
};
// 1 spanning + 4 partial = 5 reads; default min_cov = 5/2+1 = 3.
// Boundary nodes only have coverage 1 (from spanning seed) → trimmed.
let reads = vec![
spanning.clone(),
partial.clone(),
partial.clone(),
partial.clone(),
partial.clone(),
];
let result = consensus_cfg(&reads, 0, cfg);
assert!(
result.len() < 12,
"expected boundary trim with default min_cov, got len {} seq '{}'",
result.len(),
s(&result)
);
}
#[test]
fn one_spanning_many_partial_low_min_cov_reaches_partial_end() {
// Lowering min_coverage_fraction removes the boundary-trim truncation,
// but the (weight-1) normalisation in the heaviest path is a separate gate:
// edges traversed by only one read contribute score 0, and `find` prefers
// the shortest equal-score terminus. So the consensus extends to the end
// of the partial reads (position 8) but not to the spanning-only tail
// (positions 8-11) because those edges score 0.
//
// To get the full-length consensus you need ≥ 2 reads covering the tail;
// see `two_spanning_many_partial_full_length` below.
let spanning = b("ACGTACGTACGT");
let partial = b("ACGTACGT");
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let reads = vec![
spanning.clone(),
partial.clone(),
partial.clone(),
partial.clone(),
partial.clone(),
];
let result = consensus_cfg(&reads, 0, cfg);
assert!(
result.len() >= 8,
"consensus should reach at least the partial read end, got len {} seq '{}'",
result.len(),
s(&result)
);
}
#[test]
fn two_spanning_many_partial_full_length() {
// With ≥ 2 spanning reads the tail edges (nodes 8-11) get weight ≥ 2,
// contributing a positive score to the heaviest path. The boundary trim
// uses min_cov = ceil(6 * 0.1) = 1, which the spanning reads satisfy.
// Result: full 12-bp consensus.
let spanning = b("ACGTACGTACGT");
let partial = b("ACGTACGT");
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let reads = vec![
spanning.clone(),
spanning.clone(), // 2 spanning → tail edges weight=2
partial.clone(),
partial.clone(),
partial.clone(),
partial.clone(),
];
let result = consensus_cfg(&reads, 0, cfg);
assert_eq!(
result.len(),
12,
"two spanning reads should yield full-length consensus, got len {} seq '{}'",
result.len(),
s(&result)
);
}
#[test]
fn overlapping_partial_reads_assemble_beyond_seed_length() {
// Left-partial + right-partial reads that together span a longer sequence
// than any individual read. The seed covers only the left half; reads
// covering the right half extend the graph via Insert ops. With
// min_coverage_fraction = 0.1, the assembled consensus is longer than the seed.
//
// Sequence: ACGTACGTACGTACGT (16 bp)
// Left reads (12 bp): ACGTACGTACGT
// Right reads (12 bp): ACGTACGTACGT (offset 4 in the full sequence → ACGTACGTACGT)
// Together they overlap for 8 bp and cover the full 16 bp.
let left = b("ACGTACGTACGT"); // covers positions 0-11
let right = b("ACGTACGTACGT"); // same sequence; in a real scenario these
// would be from a different region, but here
// we just verify the graph can grow past seed length.
// Use an explicit short seed so the right reads extend the graph.
let seed = b("ACGTACGT"); // 8 bp seed
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
// 1 seed + 3 left (12 bp) + 3 right (12 bp) = 7 reads.
// min_cov = ceil(7 * 0.1) = 1; boundary nodes survive.
let reads = vec![
seed.clone(),
left.clone(),
left.clone(),
left.clone(),
right.clone(),
right.clone(),
right.clone(),
];
let result = consensus_cfg(&reads, 0, cfg);
assert!(
result.len() >= seed.len(),
"assembled consensus should be at least as long as the seed, got len {} seq '{}'",
result.len(),
s(&result)
);
}
/// Deterministic xorshift RNG (mirrors `tests/sv_analysis.rs`'s helper of the
/// same name) -- used only to generate a non-repetitive truth sequence below,
/// so the test isn't entangled with any of the repeat-specific known bugs.
fn xorshift_local(state: &mut u64) -> u64 {
*state ^= *state << 13;
*state ^= *state >> 7;
*state ^= *state << 17;
*state
}
#[test]
fn partial_read_population_default_coverage_floor_reaches_both_ends() {
// Regression test for the `coverage_threshold()` absolute-floor fallback
// (min_coverage_fraction == 0.0, the PoaConfig default) being sized off
// `self.n_reads` -- the *total* read count -- instead of a local
// population estimate. For a genuinely partial-read population (no
// single read spans the whole target region, by construction), that
// global floor is unreachable at either end even though each end has
// substantial, self-consistent local support, and the default-config
// consensus silently collapsed to roughly the shared middle third.
//
// Confirmed on `bench/compare_callers.py --general`'s
// `gen_short_reads_long_region_ont_r10` scenario (poa-consensus scored
// worse than both abPOA and SPOA); this is the minimal, deterministic,
// synthetic reproduction of the same population-accounting bug.
//
// Truth: 240 bp, non-repetitive (xorshift-generated, not a tandem
// repeat, to avoid entanglement with the repeat-specific known bugs).
// Left group: 8 reads covering truth[0..180] (left 180 bp).
// Right group: 8 reads covering truth[60..240] (right 180 bp).
// No read spans the full 240 bp; the two groups overlap only in the
// middle third (60..180). Before the fix: n_reads=16, old min_cov =
// 16/2+1 = 9, but each boundary region only ever has 8 reads' worth of
// agreement (the *other* group's reads never reach there at all) --
// 8 < 9, so both ends were trimmed away, leaving only the ~120 bp
// overlap. After the fix, the local population near each boundary
// correctly reflects the ~8 reads that actually reach there, and the
// full 240 bp is kept.
let mut rng: u64 = 0x5EED_1234_ABCD_0001;
let truth: Vec<u8> = (0..240)
.map(|_| b"ACGT"[(xorshift_local(&mut rng) % 4) as usize])
.collect();
let left = truth[0..180].to_vec();
let right = truth[60..240].to_vec();
assert!(left.len() < truth.len() && right.len() < truth.len());
let mut reads: Vec<Vec<u8>> = Vec::new();
for _ in 0..8 {
reads.push(left.clone());
}
for _ in 0..8 {
reads.push(right.clone());
}
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default() // min_coverage_fraction left at the 0.0 sentinel
};
let result = consensus_cfg(&reads, 0, cfg);
assert!(
result.len() >= truth.len() - 5,
"partial-read population should not collapse the default-config \
consensus to the shared middle region: got len {} vs truth len {} \
(seq: '{}')",
result.len(),
truth.len(),
s(&result)
);
// Allow a couple of bases of ordinary boundary-trim slack at the very
// tail (this test is about the population accounting, not pixel-perfect
// boundary placement); anything the fix keeps must still exactly match
// the noise-free truth base-for-base, no fabricated/rearranged content.
let n = result.len().min(truth.len());
assert_eq!(
result[..n],
truth[..n],
"recovered consensus must match the (noise-free) truth base-for-base, \
not merely in length -- got seq '{}'",
s(&result)
);
}
#[test]
fn coverage_vec_reflects_partial_read_depth() {
// The Consensus::coverage field must show lower values at positions that
// only spanning reads covered and higher values where partial reads also
// contributed.
let spanning = b("ACGTACGTACGT"); // 12 bp
let partial = b("ACGTACGT"); // 8 bp (covers prefix nodes)
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&spanning, cfg).unwrap();
for _ in 0..4 {
graph.add_read(&partial).unwrap();
}
let cons = graph.consensus().unwrap();
// First 8 positions covered by all 5 reads (spanning + 4 partial).
let prefix_min = cons.coverage[..8].iter().copied().min().unwrap_or(0);
// Last 4 positions covered only by the spanning seed.
let suffix_max = cons.coverage[8..].iter().copied().max().unwrap_or(0);
assert!(
prefix_min > suffix_max,
"prefix coverage ({}) should exceed suffix coverage ({})",
prefix_min,
suffix_max
);
}
// ── Coverage gap detection ────────────────────────────────────────────────────
#[test]
fn no_gap_when_reads_overlap() {
// Spanning seed + partials that all overlap in the middle → no coverage gap.
let seed = b("ACGTTGCAATGC"); // 12 bp
let left = b("ACGTTGCA"); // 8 bp: covers positions 0-7
let right = b("GCAATGC"); // 7 bp: covers positions 5-11 (3 bp overlap)
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..3 {
graph.add_read(&left).unwrap();
}
for _ in 0..3 {
graph.add_read(&right).unwrap();
}
let cons = graph.consensus().unwrap();
assert!(
cons.gaps.is_empty(),
"overlapping partials should produce no coverage gap; got {:?}",
cons.gaps
);
}
#[test]
fn gap_detected_when_partials_dont_overlap() {
// Spanning seed + left partials + right partials with no overlap.
// Seed: ACGTTGCAATGCCCGG (16 bp)
// Left: ACGTT (5 bp, covers positions 0-4)
// Right: CCCGG (5 bp, covers positions 11-15)
// Gap: positions 5-10 (6 bp, seed-only coverage=1)
let seed = b("ACGTTGCAATGCCCGG"); // 16 bp
let left = b("ACGTT"); // 5 bp: unique prefix of seed
let right = b("CCCGG"); // 5 bp: unique suffix of seed
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..3 {
graph.add_read(&left).unwrap();
}
for _ in 0..3 {
graph.add_read(&right).unwrap();
}
let cons = graph.consensus().unwrap();
assert!(
!cons.gaps.is_empty(),
"non-overlapping partials should produce a coverage gap; coverage={:?}",
cons.coverage
);
let gap = &cons.gaps[0];
assert!(
gap.size() >= 6,
"gap should span the 6 seed-only positions: {:?}",
gap
);
assert_eq!(gap.start + gap.size(), gap.end);
}
#[test]
fn gap_size_is_minimum_size_estimate() {
// Construct a scenario with a known gap width to verify size().
// Seed: 20 bp. Left reads cover 0-4, right reads cover 15-19.
// The middle 10 positions (5-14) have coverage=1 (seed only).
let seed = b("ACGTTGCAATGCCCGGTTAA"); // 20 bp
let left = b("ACGTT"); // 5 bp: covers 0-4
let right = b("GTTAA"); // 5 bp: covers 15-19
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..4 {
graph.add_read(&left).unwrap();
}
for _ in 0..4 {
graph.add_read(&right).unwrap();
}
let cons = graph.consensus().unwrap();
assert!(
!cons.gaps.is_empty(),
"expected a gap; coverage: {:?}",
cons.coverage
);
// The gap should span the seed-only middle region (at least 10 bp).
let total_gap: usize = cons.gaps.iter().map(|g| g.size()).sum();
assert!(
total_gap >= 10,
"expected gap ≥ 10 bp, got {total_gap}; gaps: {:?}",
cons.gaps
);
}
#[test]
fn single_read_has_no_gap() {
// With only the seed read, coverage is all 1s but there are no
// well-supported flanks, so detect_coverage_gaps returns empty.
let cfg = PoaConfig {
min_reads: 1,
..Default::default()
};
let graph = PoaGraph::new(b"ACGTTGCAATGC", cfg).unwrap();
let cons = graph.consensus().unwrap();
assert!(
cons.gaps.is_empty(),
"single-read consensus should have no gaps"
);
}
#[test]
fn gap_kind_spanning_for_seed_based_gap() {
// A seed-based gap must have kind=Spanning so callers know a minimum
// size estimate is available via size().
let seed = b("ACGTTGCAATGCCCGG");
let left = b("ACGTT");
let right = b("CCCGG");
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..3 {
graph.add_read(&left).unwrap();
}
for _ in 0..3 {
graph.add_read(&right).unwrap();
}
let cons = graph.consensus().unwrap();
assert!(!cons.gaps.is_empty());
assert_eq!(
cons.gaps[0].kind,
poa_consensus::GapKind::Spanning,
"seed-based gaps must be Spanning"
);
assert_eq!(cons.gaps[0].min_size(), Some(cons.gaps[0].size()));
}
#[test]
fn bridged_consensus_unknown_gap() {
// Two completely disjoint read groups — left reads, then right reads —
// with no read spanning the middle. bridged_consensus should produce a
// single Consensus whose gaps contain exactly one Unknown gap at the join.
let left_reads: Vec<Vec<u8>> = (0..4).map(|_| b("ACGTTGCA")).collect();
let right_reads: Vec<Vec<u8>> = (0..4).map(|_| b("ATGCCCGG")).collect();
let left_refs: Vec<&[u8]> = left_reads.iter().map(|r| r.as_slice()).collect();
let right_refs: Vec<&[u8]> = right_reads.iter().map(|r| r.as_slice()).collect();
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default()
};
let cons = poa_consensus::bridged_consensus(&left_refs, 0, &right_refs, 0, &cfg).unwrap();
// Sequence is the concatenation of both consensuses.
assert!(!cons.sequence.is_empty());
// Exactly one Unknown gap at the join point.
let unknown: Vec<_> = cons
.gaps
.iter()
.filter(|g| g.kind == poa_consensus::GapKind::Unknown)
.collect();
assert_eq!(unknown.len(), 1, "expected exactly one Unknown gap");
let gap = unknown[0];
assert_eq!(
gap.start, gap.end,
"Unknown gap should be an insertion point (start==end)"
);
assert_eq!(gap.min_size(), None, "Unknown gap has no minimum size");
// Total minimum size: at least as long as the two consensus segments.
assert!(cons.sequence.len() > 0);
assert_eq!(cons.n_reads, 8);
}
// ── path_weights and weight_fraction ─────────────────────────────────────────
#[test]
fn path_weights_reflect_edge_support() {
// 1 spanning seed + 4 partial reads covering the first 8 of 12 nodes.
// Interior edges (0-7) should have weight 5 (seed + 4 partial).
// The partial reads end at node 7, so the consensus (heaviest path) stops
// there. All 8 weights should be ≥ 2 (shared) and n_reads should be 5.
let spanning = b("ACGTACGTACGT");
let partial = b("ACGTACGT");
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&spanning, cfg).unwrap();
for _ in 0..4 {
graph.add_read(&partial).unwrap();
}
let cons = graph.consensus().unwrap();
assert_eq!(cons.n_reads, 5);
assert_eq!(cons.path_weights.len(), cons.sequence.len());
// Every weight in the consensus should reflect multi-read support.
for (i, &w) in cons.path_weights.iter().enumerate() {
assert!(
w >= 2,
"position {i}: weight {w} should be ≥ 2 (shared by seed + partial)"
);
}
}
#[test]
fn weight_fraction_in_unit_interval() {
let reads = vec![b("ACGTACGT"); 5];
let cons = consensus_cfg(&reads, 0, PoaConfig::default());
// Build Consensus directly to check the fraction helper.
let mut graph = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
for r in &reads[1..] {
graph.add_read(r).unwrap();
}
let c = graph.consensus().unwrap();
let fracs = c.weight_fraction();
assert_eq!(fracs.len(), cons.len());
for (i, &f) in fracs.iter().enumerate() {
assert!(
(0.0..=1.0).contains(&f),
"position {i}: fraction {f} out of [0,1]"
);
}
// All reads identical → all fractions should be 1.0.
for (i, &f) in fracs.iter().enumerate() {
assert!(
(f - 1.0).abs() < 1e-6,
"position {i}: expected fraction 1.0, got {f}"
);
}
}
#[test]
fn weight_fraction_drops_at_single_read_positions() {
// Non-repetitive spanning sequence so partial reads have exactly one valid
// alignment position (avoids the rotation-phase tie-break issue that arises
// with periodic sequences like ACGTACGTACGT).
//
// 2 spanning reads (needed so tail edges score > 0 and enter the path) +
// 4 partial reads covering only the first 8 of 12 bases.
// Tail positions (8-11) are supported only by the 2 spanning reads;
// their fraction (2/6 ≈ 0.33) should be below the prefix fraction (6/6 = 1.0).
let spanning = b("ACGTTGCAATGC"); // 12 bp, no 8-mer repeats
let partial = b("ACGTTGCA"); // 8 bp, uniquely matches positions 0-7
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
min_coverage_fraction: 0.1,
..Default::default()
};
let mut graph = PoaGraph::new(&spanning, cfg.clone()).unwrap();
graph.add_read(&spanning).unwrap();
for _ in 0..4 {
graph.add_read(&partial).unwrap();
}
let cons = graph.consensus().unwrap();
let fracs = cons.weight_fraction();
assert_eq!(cons.sequence.len(), 12, "expected full-length consensus");
let prefix_frac: f32 = fracs[..8].iter().copied().sum::<f32>() / 8.0;
let suffix_frac: f32 = fracs[8..].iter().copied().sum::<f32>() / 4.0;
assert!(
prefix_frac > suffix_frac,
"prefix avg fraction ({prefix_frac:.2}) should exceed suffix ({suffix_frac:.2})"
);
}
// Semi-global op-level tests: verify the alignment ops themselves, not just
// the consensus. These catch bugs that happen to not affect the output length
// but still corrupt edge weights or delete_counts.
#[test]
fn semi_global_no_prefix_deletes_for_mid_start_read() {
// Seed "GGACGT", partial read "ACGT" matches the suffix perfectly.
// In global mode the aligner is forced to start from the source (G,G) and
// emits Delete ops for those prefix nodes. Semi-global lets the read start
// at the first matching node and must produce zero Deletes.
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
band_width: 0,
..Default::default()
};
let graph = PoaGraph::new(b"GGACGT", cfg).unwrap();
let (ops, _) = graph.align_read_ops_unbanded(b"ACGT").unwrap();
let n_del = ops
.iter()
.filter(|op| matches!(op, AlignOp::Delete(_)))
.count();
let n_mat = ops
.iter()
.filter(|op| matches!(op, AlignOp::Match(_)))
.count();
assert_eq!(
n_del, 0,
"semi-global: expected no prefix Deletes, got {:?}",
ops
);
assert_eq!(n_mat, 4, "semi-global: expected 4 Matches, got {:?}", ops);
}
#[test]
fn global_produces_prefix_deletes_for_mid_start_read() {
// Same setup, global mode: the alignment is forced through the GG prefix,
// producing Delete ops for those nodes.
let cfg = PoaConfig {
band_width: 0,
alignment_mode: AlignmentMode::Global,
..Default::default()
};
let graph = PoaGraph::new(b"GGACGT", cfg).unwrap();
let (ops, _) = graph.align_read_ops_unbanded(b"ACGT").unwrap();
let n_del = ops
.iter()
.filter(|op| matches!(op, AlignOp::Delete(_)))
.count();
assert!(
n_del > 0,
"global: expected prefix Delete ops for mid-start read"
);
}
#[test]
fn semi_global_spanning_read_matches_global() {
// A read that spans the full seed is not partial; semi-global and global
// must produce the same consensus.
let reads = vec![b("ACGTACGT"); 4];
let global = consensus(&reads, 0);
let cfg = PoaConfig {
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default()
};
let semi = consensus_cfg(&reads, 0, cfg);
assert_eq!(
global, semi,
"spanning reads: semi-global must equal global"
);
}
// ── heaviest_path: Match/Delete edge-weight conflation ───────────────────────
//
// `Node.out_edges` used to store one conflated `i32` weight incremented
// identically by Match and Delete traversals (see `increment_or_add_edge`),
// and `heaviest_path`'s cumulative-weight DP summed that raw total. A node
// reached mostly by Delete (reads that skip past it, confirming nothing
// about its base) could therefore out-compete a genuinely Match-confirmed
// alternative arm at the same fork, purely on raw traffic. Confirmed on real
// RFC1 CANVAS data during the architectural audit (see
// design/graph_data_model_rework.md); this test is the minimal, permanent,
// forced reproduction of that same mechanism.
#[test]
fn heaviest_path_prefers_matched_over_delete_inflated_arm() {
// Non-repetitive flanks: semi-global's free leading/trailing gap can't
// "explain away" the length deficit of the deletion allele by shifting
// a homopolymer-adjacent frame for free (that alternative alignment
// would incur many mismatches against this flank), so the aligner is
// forced to represent the deletion allele's reads as a genuine interior
// Delete rather than a boundary-shift mismatch.
let left = b"ACGTGCAT";
let right = b"TACGATCG";
let mut seed = Vec::new();
seed.extend_from_slice(left);
seed.push(b'G'); // "reference" base at the contested position
seed.extend_from_slice(right);
let mut del_read = Vec::new(); // true deletion allele: missing the middle base
del_read.extend_from_slice(left);
del_read.extend_from_slice(right);
let mut snp_read = Vec::new(); // true SNP allele: G -> C at the same position
snp_read.extend_from_slice(left);
snp_read.push(b'C');
snp_read.extend_from_slice(right);
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();
// 6 reads: true deletion allele -- should register as genuine interior
// Delete ops against the seed's 'G' node, not as Match.
for _ in 0..6 {
g.add_read(&del_read).unwrap();
}
// 4 reads: true SNP allele -- genuine Match confirmation of a 'C' node.
for _ in 0..4 {
g.add_read(&snp_read).unwrap();
}
let topo = g.graph_topology();
let g_node = topo
.nodes
.iter()
.find(|n| n.base == b'G' && n.coverage == 1 && n.delete_count == 6)
.expect("expected the reference 'G' node with cov=1, delete_count=6 (6 genuine interior deletes)");
let c_node = topo
.nodes
.iter()
.find(|n| n.base == b'C' && n.coverage == 4 && n.delete_count == 0)
.expect("expected the SNP 'C' node with cov=4, delete_count=0 (4 genuine matches)");
// Sanity: confirm the setup itself is what it claims to be before
// asserting anything about which one heaviest_path picked. If this
// fails, the *test* is wrong (e.g. the aligner took a different path
// than the forced-delete design assumes), not the fix.
assert_eq!(
g_node.coverage, 1,
"'G' should have only the seed's own match"
);
assert_eq!(
g_node.delete_count, 6,
"'G' should be deleted-through by all 6 deletion-allele reads"
);
assert_eq!(
c_node.coverage, 4,
"'C' should be matched by all 4 SNP-allele reads"
);
assert_eq!(c_node.delete_count, 0, "'C' is never deleted through");
let spine_ranks: std::collections::HashSet<usize> = topo.spine_ranks.iter().copied().collect();
assert!(
spine_ranks.contains(&c_node.topo_rank),
"heaviest_path should select the genuinely Match-confirmed 'C' node (4 real \
confirmations) over the Delete-inflated 'G' node (1 real confirmation, 6 skips) -- \
it did not, meaning delete traffic is still winning the routing decision"
);
assert!(
!spine_ranks.contains(&g_node.topo_rank),
"'G' -- confirmed by only the seed itself, with 6 reads skipping past it -- \
should not be on the heaviest-path spine ahead of the better-evidenced 'C'"
);
}
// ── edge_reads: Match/Delete traversal-type split (Phase 2) ─────────────────
//
// `PoaGraph.edge_reads` used to record read membership for ANY traversal type
// (Match, founding Insert, or Delete) with no distinction -- so a read that
// merely deleted through a bubble arm's starting node (skipped it, confirming
// nothing) was indistinguishable from a read that genuinely matched that arm.
// `partition_reads_by_bubble`, `phasing_groups`, and `BubbleSite.arm_read_counts`
// all read `edge_reads` to decide "which reads support which arm" -- so all
// three inherited the same blindness. Delete traversals now go into a separate
// `edge_delete_reads` map instead (see design/graph_data_model_rework.md Phase 2).
//
// Both tests below share the same forced setup: non-repetitive flanks (so
// semi-global's free boundary gap can't explain the length deficit for free,
// forcing a genuine interior Delete -- see the Phase 1 regression test for the
// same technique), a SNP-style 1-node bubble (G vs C) with the two arms' raw
// read counts deliberately near-tied (3 `g_read`s + the seed's own base = 4
// total for G; 4 `c_read`s = 4 total for C) before adding one more read that
// is missing the contested base entirely. Read indices: seed=0, g_reads=1-3,
// c_reads=4-7, the deletion read=8.
fn phase2_bubble_setup() -> (Vec<u8>, Vec<u8>, Vec<u8>, Vec<u8>) {
let left = b"ACGTGCAT";
let right = b"TACGATCG";
let mut seed = Vec::new();
seed.extend_from_slice(left);
seed.push(b'G'); // "reference" base at the contested position
seed.extend_from_slice(right);
let mut g_read = Vec::new();
g_read.extend_from_slice(left);
g_read.push(b'G');
g_read.extend_from_slice(right);
let mut c_read = Vec::new();
c_read.extend_from_slice(left);
c_read.push(b'C'); // true SNP allele
c_read.extend_from_slice(right);
let mut del_read = Vec::new(); // missing the contested base entirely
del_read.extend_from_slice(left);
del_read.extend_from_slice(right);
(seed, g_read, c_read, del_read)
}
fn phase2_bubble_graph() -> PoaGraph {
let (seed, g_read, c_read, del_read) = phase2_bubble_setup();
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2,
band_width: 0,
adaptive_band: false,
warn_on_long_unbanded: false,
..Default::default()
};
let mut g = PoaGraph::new(&seed, cfg).unwrap();
for _ in 0..3 {
g.add_read(&g_read).unwrap();
}
for _ in 0..4 {
g.add_read(&c_read).unwrap();
}
g.add_read(&del_read).unwrap();
g
}
/// `partition_reads_by_bubble` / `phasing_groups` (exercised here via
/// `consensus_multi`) must not attribute a read to an arm it merely deleted
/// through -- only reads that genuinely confirmed an arm's base should
/// determine group membership.
#[test]
fn partition_reads_by_bubble_excludes_delete_only_reads() {
let g = phase2_bubble_graph();
// Sanity: confirm the setup produced the fork shape this test assumes,
// before asserting anything about the fix. If this fails, the *test* is
// wrong (e.g. the aligner didn't force a genuine interior Delete the way
// the non-repetitive-flank design assumes), not the fix.
let topo = g.graph_topology();
let g_node = topo
.nodes
.iter()
.find(|n| n.base == b'G' && n.coverage == 4 && n.delete_count == 1)
.expect(
"expected 'G' with coverage=4 (seed + 3 g_reads), delete_count=1 (the deletion read)",
);
let c_node = topo
.nodes
.iter()
.find(|n| n.base == b'C' && n.coverage == 4 && n.delete_count == 0)
.expect("expected 'C' with coverage=4 (4 c_reads), delete_count=0");
assert_eq!(g_node.coverage, 4);
assert_eq!(c_node.coverage, 4);
let multi = g.consensus_multi().unwrap();
assert_eq!(multi.len(), 2, "expected two allele consensuses (G and C)");
let g_allele = multi
.iter()
.find(|c| c.sequence.contains(&b'G') && c.read_indices.contains(&1))
.expect("expected to find the G-allele consensus (contains read 1, a g_read)");
let c_allele = multi
.iter()
.find(|c| c.read_indices.contains(&4))
.expect("expected to find the C-allele consensus (contains read 4, a c_read)");
assert!(
!g_allele.read_indices.contains(&8),
"read 8 deleted through 'G' -- it must not be attributed to the G arm's group just \
because it structurally traversed the same edge; got read_indices={:?}",
g_allele.read_indices
);
// Read 8 has no genuine confirmation of either arm, so it is legitimately
// "unassigned" and folds into whichever group is largest (C, at 4 reads
// vs G's 3) -- that's a defensible default for an ambiguous read, not a
// wrong attribution the way counting it as a G-confirmation would be.
assert!(
c_allele.read_indices.contains(&8),
"expected the unassigned deletion read (8) to fold into the larger (C) group; \
got read_indices={:?}",
c_allele.read_indices
);
assert_eq!(
g_allele.read_indices,
vec![0, 1, 2, 3],
"G's group should be exactly the seed + 3 g_reads, no more"
);
}
/// `BubbleSite.arm_read_counts` must reflect only genuine Match-confirmed
/// reads for each arm, not reads that merely deleted through its entry node.
#[test]
fn bubble_site_arm_read_counts_excludes_delete_only_reads() {
let g = phase2_bubble_graph();
let cons = g.consensus().unwrap();
assert_eq!(
cons.bubble_sites.len(),
1,
"expected exactly one bubble site"
);
let site = &cons.bubble_sites[0];
assert_eq!(site.arm_sequences.len(), 2);
let g_idx = site
.arm_sequences
.iter()
.position(|s| s == b"G")
.expect("expected a 'G' arm");
let c_idx = site
.arm_sequences
.iter()
.position(|s| s == b"C")
.expect("expected a 'C' arm");
assert_eq!(
site.arm_read_counts[g_idx], 4,
"'G' arm should count only its 4 genuine confirmations (seed + 3 g_reads), \
not the 1 additional read that merely deleted through it (would be 5 if \
still conflated); got {:?}",
site.arm_read_counts
);
assert_eq!(
site.arm_read_counts[c_idx], 4,
"'C' arm's count is unaffected either way (no deletes through it); got {:?}",
site.arm_read_counts
);
}
// ── Arm content-addressing (Phase 3) ─────────────────────────────────────────
//
// Re-run of the period-7 (`GCTAGCT`x10) duplicate-fork audit from the
// architectural investigation that motivated design/graph_data_model_rework.md.
// A tandem-repeat context produces many small, per-read insertion/substitution
// arms; before Phase 3, exact-duplicate arms at the same fork (same base(s),
// same position) still fragmented into separate `coverage=1` singleton nodes
// unless ordinary DP happened to notice the match on its own. Content
// addressing (a per-fork `HashMap<Vec<u8>, node_idx>` keyed by the *full*
// characterized edit, checked before creating a new node) collapses genuine
// exact duplicates into one shared, higher-coverage node.
//
// Measured directly (apples-to-apples: identical simulated read set, same
// unbanded config, compared against an isolated build of the Phase 1+2 commit
// with no Phase 3 code at all):
// before (Phase 1+2 only): 166 nodes, 71 singleton(cov=1) nodes, 14 duplicate-fork positions
// after (Phase 3): 163 nodes, 67 singleton(cov=1) nodes, 12 duplicate-fork positions
//
// This is a real but modest reduction, not a collapse -- consistent with the
// design doc's own prediction that this specific scenario is dominated by
// *fuzzy* near-duplicate arms (same effective edit, different incidental
// length/shape, e.g. a 4-base vs. a 5-base insertion reconverging at the same
// point) rather than byte-identical exact duplicates. Fuzzy near-duplicates
// are explicitly out of scope for Phase 3 (would need canonicalization, a
// follow-on problem per the design doc) and are NOT expected to merge here.
//
// **Pure-bypass (bypass_edge_delete_rework revised Phase 1) further reduced
// the structural counts to 162 nodes / 65 singletons / 12 duplicate forks.**
// This is a *representation* change, not an output change: the deleting reads
// no longer create delete-bucket edges into skipped nodes, and post-delete
// divergences hang off the entry predecessor rather than the last deleted
// node, so a couple of fragmentation artifacts disappear -- in the same
// simplifying direction as the whole rework. The consensus OUTPUT is
// unchanged at 77bp (one extra repeat unit over the 70bp truth) across all of
// Phase 1+2, Phase 3, AND pure bypass -- asserted explicitly below as the
// actual correctness invariant, so the structural counts are documented drift,
// not a magic-number gate on behavior.
#[test]
fn period7_content_addressing_reduces_duplicate_forks() {
let template = repeat_unit(b"GCTAGCT", 10); // 70 bp
let reads = simulate_reads(&template, 20, 0.05, 0.02, 0.02, 18);
let cfg = PoaConfig {
band_width: 0,
adaptive_band: false,
warn_on_long_unbanded: false,
..Default::default()
};
let mut g = PoaGraph::new(&reads[0], cfg).unwrap();
for read in &reads[1..] {
g.add_read(read).unwrap();
}
let topo = g.graph_topology();
use std::collections::HashMap;
let base_of: HashMap<usize, u8> = topo.nodes.iter().map(|n| (n.topo_rank, n.base)).collect();
let cov_of: HashMap<usize, u32> = topo
.nodes
.iter()
.map(|n| (n.topo_rank, n.coverage))
.collect();
let mut by_from: HashMap<usize, Vec<(u8, u32)>> = HashMap::new();
for e in &topo.edges {
by_from
.entry(e.from_rank)
.or_default()
.push((base_of[&e.to_rank], cov_of[&e.to_rank]));
}
let mut dup_fork_positions = 0usize;
for succs in by_from.values() {
if succs.len() < 2 {
continue;
}
let mut by_base: HashMap<u8, usize> = HashMap::new();
for &(b, _) in succs {
*by_base.entry(b).or_default() += 1;
}
dup_fork_positions += by_base.values().filter(|&&count| count >= 2).count();
}
let singletons = topo.nodes.iter().filter(|n| n.coverage == 1).count();
// Output-correctness invariant (unchanged across Phase 1+2, Phase 3, and
// pure bypass): the consensus is 77bp, one extra repeat unit over the 70bp
// truth. This is the property that must NOT move; the structural counts
// below are documented representation drift.
let cons_len = g.consensus().unwrap().sequence.len();
assert_eq!(
cons_len, 77,
"consensus output length must stay 77bp (the scenario's invariant); \
a change here is an output-correctness regression, not representation drift"
);
// Structural fragmentation metrics -- representation, not output. Updated
// to the pure-bypass baseline (was 163/67/12 at Phase 3; see the module
// comment above for why each dropped).
assert_eq!(
topo.nodes.len(),
162,
"node count drifted from the pure-bypass baseline (163 at Phase 3, 162 under pure bypass)"
);
assert_eq!(
singletons, 65,
"singleton(cov=1) node count drifted from the pure-bypass baseline (67 at Phase 3, 65 under pure bypass)"
);
assert_eq!(
dup_fork_positions, 12,
"duplicate-fork-position count drifted from the pure-bypass baseline (12 at Phase 3, 12 under pure bypass)"
);
}
fn xorshift(state: &mut u64) -> u64 {
*state ^= *state << 13;
*state ^= *state >> 7;
*state ^= *state << 17;
*state
}
fn rand_f64(state: &mut u64) -> f64 {
xorshift(state) as f64 / u64::MAX as f64
}
fn random_base(state: &mut u64) -> u8 {
b"ACGT"[(xorshift(state) % 4) as usize]
}
fn random_base_not(exclude: u8, state: &mut u64) -> u8 {
let opts: [u8; 3] = match exclude {
b'A' => [b'C', b'G', b'T'],
b'C' => [b'A', b'G', b'T'],
b'G' => [b'A', b'C', b'T'],
_ => [b'A', b'C', b'G'],
};
opts[(xorshift(state) % 3) as usize]
}
fn simulate_reads(
template: &[u8],
n_reads: usize,
sub: f64,
ins: f64,
del: f64,
seed: u64,
) -> Vec<Vec<u8>> {
let mut state = seed;
let mut reads = Vec::with_capacity(n_reads);
for _ in 0..n_reads {
let mut read = Vec::with_capacity(template.len());
for &base in template {
if rand_f64(&mut state) < ins {
read.push(random_base(&mut state));
}
if rand_f64(&mut state) < del {
continue;
}
if rand_f64(&mut state) < sub {
read.push(random_base_not(base, &mut state));
} else {
read.push(base);
}
}
if !read.is_empty() {
reads.push(read);
}
}
reads
}
fn repeat_unit(unit: &[u8], n: usize) -> Vec<u8> {
unit.iter().cycle().take(unit.len() * n).copied().collect()
}
// ── Reverse complement / orientation ─────────────────────────────────────────
#[test]
fn reverse_complement_basic() {
use crate::reverse_complement;
assert_eq!(reverse_complement(b"ACGT"), b"ACGT");
assert_eq!(reverse_complement(b"AAAA"), b"TTTT");
assert_eq!(reverse_complement(b"GCTA"), b"TAGC");
}
#[test]
fn orient_to_seed_forward() {
use crate::Strand;
use crate::orient_to_seed;
let seed = b("ACGTACGTACGT");
let read = b("ACGTACGT");
assert_eq!(orient_to_seed(&read, &seed, 4), Strand::Forward);
}
#[test]
fn orient_to_seed_reverse() {
use crate::Strand;
use crate::orient_to_seed;
// Use a non-palindromic sequence so forward and RC share no k-mers.
let seed = b("AAAACCCCGGGG");
let rc = crate::reverse_complement(&seed);
assert_eq!(orient_to_seed(&rc, &seed, 4), Strand::Reverse);
}
#[test]
fn mixed_strand_input() {
use crate::auto_orient;
let seed = b("CATCATCAT");
let rc = crate::reverse_complement(&seed);
let reads = vec![seed.clone(), seed.clone(), rc.clone(), rc.clone()];
let oriented: Vec<Vec<u8>> = auto_orient(&reads, 0)
.into_iter()
.map(|c| c.into_owned())
.collect();
// All oriented reads should match or be rc'd to match seed strand.
let mut graph = PoaGraph::new(&oriented[0], PoaConfig::default()).unwrap();
for r in &oriented[1..] {
graph.add_read(r).unwrap();
}
let result = graph.consensus().unwrap().sequence;
assert_eq!(result.len(), 9, "mixed strand: got len {}", result.len());
}
// ── Majority-frequency consensus ──────────────────────────────────────────────
fn mf_cfg() -> PoaConfig {
PoaConfig {
consensus_mode: ConsensusMode::MajorityFrequency,
..Default::default()
}
}
fn mf_consensus(reads: &[Vec<u8>], seed_idx: usize) -> Vec<u8> {
consensus_cfg(reads, seed_idx, mf_cfg())
}
#[test]
fn mf_identical_reads() {
let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCAT")];
assert_eq!(mf_consensus(&reads, 0), b("CATCATCAT"));
}
#[test]
fn mf_matches_hb_on_clean_input() {
// On a read set with no noise, MF and HB should agree.
let reads = vec![
b("CAGCAGCAG"),
b("CAGCAGCAG"),
b("CAGCAGCAGCAG"),
b("CAGCAGCAG"),
];
let hb = consensus(&reads, 0);
let mf = mf_consensus(&reads, 0);
assert_eq!(hb, mf, "HB and MF disagree on clean input");
}
#[test]
fn mf_boundary_trim_leading() {
// Seed has 3 extra leading bases not present in the majority.
// MF should exclude them because gap votes outnumber base votes.
let reads = vec![
b("XXXCATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let result = s(&mf_consensus(&reads, 0));
assert_eq!(result, "CATCATCAT", "got: {}", result);
}
#[test]
fn mf_boundary_trim_trailing() {
let reads = vec![
b("CATCATCATXXX"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let result = s(&mf_consensus(&reads, 0));
assert_eq!(result, "CATCATCAT", "got: {}", result);
}
#[test]
fn mf_majority_base_wins() {
// 3 reads have CAT, 1 has CGT at position 1. MF should pick A.
let reads = vec![
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CGTCATCAT"),
];
assert_eq!(s(&mf_consensus(&reads, 0)), "CATCATCAT");
}
#[test]
fn mf_single_outlier_not_inflated() {
// One read has an extra CAT; MF should not include it.
let reads = vec![b("CATCATCAT"), b("CATCATCAT"), b("CATCATCATCAT")];
assert_eq!(mf_consensus(&reads, 0).len(), 9);
}
// ── GraphStats ────────────────────────────────────────────────────────────────
fn build_graph(reads: &[Vec<u8>], seed_idx: usize) -> PoaGraph {
let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
for (i, read) in reads.iter().enumerate() {
if i != seed_idx {
graph.add_read(read).unwrap();
}
}
graph
}
#[test]
fn stats_clean_linear_no_bubbles() {
// Identical reads produce a clean linear graph with no bubbles.
let reads = vec![
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let st = build_graph(&reads, 0).stats();
assert_eq!(st.bubble_count, 0);
assert_eq!(st.max_bubble_depth, 0);
assert_eq!(st.node_count, 9);
// All reads match every node: delete_count=0 everywhere → entropy=0.
assert_eq!(st.mean_column_entropy, 0.0);
}
#[test]
fn stats_bubble_detected() {
// 3 reads with CATCATCAT, 1 with CGTCATCAT → SNV bubble at position 1 (A vs G).
let reads = vec![
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CGTCATCAT"),
];
let st = build_graph(&reads, 0).stats();
assert_eq!(st.bubble_count, 1, "expected 1 bubble");
// Minority arm has 1 read (the CGT read creates a G branch at position 1).
assert_eq!(st.max_bubble_depth, 1, "minority arm weight should be 1");
}
#[test]
fn stats_entropy_nonzero_on_length_variation() {
// Seed has 3 leading X nodes that other reads delete.
// The X nodes get delete_count=3, coverage=1 → entropy > 0.
// (Shorter reads aligned to a longer graph in global mode don't generate trailing
// deletes — they simply stop at the best-scoring diagonal. Leading deletes DO fire
// because the traceback reaches t=0, j=0 via the D-chain at the j=0 column.)
// Global alignment is required here: with SemiGlobal, shorter reads take a free
// terminal gap over the XXX prefix and never traverse those nodes, so delete_count=0.
let reads = vec![
b("XXXCATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let cfg = PoaConfig {
alignment_mode: AlignmentMode::Global,
band_width: 0,
adaptive_band: false,
..PoaConfig::default()
};
let mut graph = PoaGraph::new(&reads[0], cfg).unwrap();
for read in &reads[1..] {
graph.add_read(read).unwrap();
}
let st = graph.stats();
// X nodes: coverage=1, delete_count=3 → p=0.25, binary_entropy(0.25) ≈ 0.811 bits.
assert!(
st.mean_column_entropy > 0.0,
"expected nonzero entropy, got {}",
st.mean_column_entropy
);
}
#[test]
fn stats_node_edge_counts() {
// 2 reads with length variation: 9-node backbone + 3 extra nodes from longer read.
// The longer read aligns as Insert(CAT) + Match(nodes 0-8), creating nodes 9,10,11
// with edges 9→10, 10→11, 11→0. Backbone edges 0→1..7→8 already existed = 8.
let reads = vec![b("CATCATCAT"), b("CATCATCATCAT")];
let st = build_graph(&reads, 0).stats();
assert_eq!(st.node_count, 12, "9 + 3 extra nodes");
// 8 backbone + 3 new edges (9→10, 10→11, 11→0) = 11 total.
assert_eq!(st.edge_count, 11);
}
#[test]
fn stats_coverage_mean_uniform() {
// All 4 reads match all 9 nodes → coverage=4 everywhere → variance=0.
let reads = vec![
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
b("CATCATCAT"),
];
let st = build_graph(&reads, 0).stats();
assert!((st.coverage_mean - 4.0).abs() < 1e-10);
assert!(st.coverage_variance < 1e-10);
}
// ── Multi-allele consensus ────────────────────────────────────────────────────
fn multi_graph(reads: &[Vec<u8>], seed_idx: usize) -> PoaGraph {
let mut graph = PoaGraph::new(&reads[seed_idx], PoaConfig::default()).unwrap();
for (i, read) in reads.iter().enumerate() {
if i != seed_idx {
graph.add_read(read).unwrap();
}
}
graph
}
#[test]
fn consensus_multi_single_allele() {
// Identical reads → no bubble → consensus_multi falls through to single consensus.
let reads = vec![b("CATCATCAT"); 4];
let g = multi_graph(&reads, 0);
let results = g.consensus_multi().unwrap();
assert_eq!(results.len(), 1, "expected 1 allele for homozygous input");
assert_eq!(results[0].sequence, b("CATCATCAT"));
}
#[test]
fn consensus_multi_snv_bubble() {
// 4 reads with allele A (CATCATCAT) and 4 with allele B (CATCGTCAT).
// A SNV at position 4 (A→G) creates a 2-arm bubble.
let allele_a = b("CATCATCAT");
let allele_b = b("CATCGTCAT");
let reads: Vec<Vec<u8>> = (0..4)
.map(|_| allele_a.clone())
.chain((0..4).map(|_| allele_b.clone()))
.collect();
let g = multi_graph(&reads, 0);
let results = g.consensus_multi().unwrap();
assert_eq!(results.len(), 2, "expected 2 alleles for SNV input");
let seqs: Vec<String> = results.iter().map(|c| s(&c.sequence)).collect();
assert!(
seqs.iter().any(|seq| seq == "CATCATCAT"),
"missing CATCATCAT allele: {:?}",
seqs
);
assert!(
seqs.iter().any(|seq| seq == "CATCGTCAT"),
"missing CATCGTCAT allele: {:?}",
seqs
);
}
#[test]
fn consensus_multi_length_variation() {
// Two alleles with different repeat counts flanked by matching anchors.
// The anchor regions create a proper bubble between the two allele lengths.
// short: AAA + 2×CAT + TTTTTT = 14 bp
// long : AAA + 3×CAT + TTTTTT = 17 bp
let short = b("AAACATCATTTTTT");
let long_ = b("AAACATCATCATTTTTT");
let reads: Vec<Vec<u8>> = (0..4)
.map(|_| short.clone())
.chain((0..4).map(|_| long_.clone()))
.collect();
let g = multi_graph(&reads, 0);
let results = g.consensus_multi().unwrap();
assert_eq!(
results.len(),
2,
"expected 2 alleles for length-variation input"
);
let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
assert!(lens.contains(&14), "expected 14-bp allele; got {:?}", lens);
assert!(lens.contains(&17), "expected 17-bp allele; got {:?}", lens);
}
#[test]
fn consensus_multi_insufficient_depth_per_allele() {
// 4 reads total, min_reads=3. Each allele only gets 2 → InsufficientDepth.
let allele_a = b("CATCATCAT");
let allele_b = b("CATCGTCAT");
let reads = vec![allele_a.clone(), allele_a, allele_b.clone(), allele_b];
let cfg = PoaConfig {
min_reads: 3,
..Default::default()
};
let mut g = PoaGraph::new(&reads[0], cfg).unwrap();
for r in &reads[1..] {
g.add_read(r).unwrap();
}
let result = g.consensus_multi();
assert!(
matches!(result, Err(PoaError::InsufficientDepth { .. })),
"expected InsufficientDepth, got {:?}",
result.map(|v| v.len())
);
}
// ── Longer-sequence stress tests ──────────────────────────────────────────────
#[test]
fn long_repeat_consensus_correctness() {
let seq: Vec<u8> = "CAT".repeat(30).into_bytes(); // 90 bp
let reads = vec![seq.clone(); 6];
assert_eq!(consensus(&reads, 0), seq, "30×CAT consensus mismatch");
}
#[test]
fn long_repeat_length_majority_wins() {
// 8 reads at 60 bp (20×CAT), 2 outliers at 63 bp (21×CAT)
let maj: Vec<u8> = "CAT".repeat(20).into_bytes();
let out: Vec<u8> = "CAT".repeat(21).into_bytes();
let mut reads: Vec<Vec<u8>> = vec![maj.clone(); 8];
reads.extend(vec![out; 2]);
let result = consensus(&reads, 0);
assert_eq!(
result.len(),
60,
"expected 60-bp majority, got {} bp",
result.len()
);
}
#[test]
fn long_repeat_snv_correction() {
// 9 correct reads + 1 noisy read with a single mismatch at position 30
let correct: Vec<u8> = "CAT".repeat(20).into_bytes(); // 60 bp
let mut noisy = correct.clone();
noisy[30] = b'G';
let mut reads: Vec<Vec<u8>> = vec![correct.clone(); 9];
reads.push(noisy);
let result = consensus(&reads, 0);
assert_eq!(
result, correct,
"SNV from single noisy read should not affect consensus"
);
}
#[test]
fn long_banded_matches_unbanded() {
// 4 × 72 bp + 1 × 78 bp (length outlier), band=30
let base: Vec<u8> = "CAT".repeat(24).into_bytes(); // 72 bp
let long: Vec<u8> = "CAT".repeat(26).into_bytes(); // 78 bp
let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
reads.push(long);
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 30,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded, unbanded,
"banded(30) should match unbanded for small divergence"
);
}
#[test]
fn long_adaptive_band_matches_unbanded() {
// 4 × 72 bp + 1 × 81 bp, adaptive band b=10 f=0.05
let base: Vec<u8> = "CAT".repeat(24).into_bytes(); // 72 bp
let long: Vec<u8> = "CAT".repeat(27).into_bytes(); // 81 bp
let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
reads.push(long);
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
adaptive_band: true,
adaptive_band_b: 10,
adaptive_band_f: 0.05,
..Default::default()
};
let adaptive = consensus_cfg(&reads, 0, cfg);
assert_eq!(
adaptive, unbanded,
"adaptive band should match unbanded for small divergence"
);
}
// ── Performance optimisation tests ───────────────────────────────────────────
#[test]
fn skip_fires_on_clean_reads() {
// Zero-error reads: the diagonal skip fires on ~100% of rows.
// Verify correctness — narrow band with skip active must produce the same
// consensus as unbanded and must equal the read itself.
let read: Vec<u8> = "CAT".repeat(10).into_bytes(); // 30 bp, 10 identical reads
let reads: Vec<Vec<u8>> = vec![read.clone(); 10];
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 5,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded, unbanded,
"diagonal skip: banded must match unbanded on identical reads"
);
assert_eq!(
banded, read,
"diagonal skip: consensus of identical reads must equal the read"
);
}
#[test]
fn tracking_band_survives_phase_shift() {
// One read with a 10-bp prefix insertion shifts the alignment diagonal by 10.
// With a fixed-diagonal band of 5 this would be BandTooNarrow; with the
// tracking band the band re-centres after the shift, and smart retry widens
// the initial band so alignment succeeds.
let base: Vec<u8> = "ACGT".repeat(15).into_bytes(); // 60 bp
let shifted: Vec<u8> = {
let mut s = b"AAAAAAAAAA".to_vec(); // 10 bp prefix → diagonal drift +10
s.extend_from_slice(&base);
s
}; // 70 bp
let mut reads: Vec<Vec<u8>> = vec![base.clone(); 4];
reads.push(shifted);
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 5,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded, unbanded,
"tracking band: must match unbanded on reads with a large phase shift"
);
}
#[test]
fn sv_retry_correct() {
// One outlier read with a 15-bp expansion against 5 short reads. band_width=3
// is too narrow to track the diagonal shift: approaching-edge fires on every
// row (right_margin = w = 3 < GAP_MARGIN), forcing a smart retry with a wider
// band (~17). The retry band covers j=l=24 at the last graph node, so the
// alignment succeeds without a second retry. Majority consensus is "CAT"×3.
let short: Vec<u8> = "CAT".repeat(3).into_bytes(); // 9 bp
let expanded: Vec<u8> = "CAT".repeat(8).into_bytes(); // 24 bp (+15 bp expansion)
let mut reads: Vec<Vec<u8>> = vec![short.clone(); 5];
reads.push(expanded);
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 3,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded, unbanded,
"sv_retry: smart retry must produce correct consensus when SV read forces band widening"
);
}
#[test]
fn consensus_multi_long_flanked_str() {
// Two alleles anchored by GGGGG / AAAAA flanks:
// allele_a: GGGGG + 8×CAT + AAAAA = 5+24+5 = 34 bp
// allele_b: GGGGG + 11×CAT + AAAAA = 5+33+5 = 43 bp
let flank_l = b("GGGGG");
let flank_r = b("AAAAA");
let inner_a: Vec<u8> = "CAT".repeat(8).into_bytes();
let inner_b: Vec<u8> = "CAT".repeat(11).into_bytes();
let allele_a: Vec<u8> = [flank_l.as_slice(), inner_a.as_slice(), flank_r.as_slice()].concat();
let allele_b: Vec<u8> = [flank_l.as_slice(), inner_b.as_slice(), flank_r.as_slice()].concat();
let mut reads: Vec<Vec<u8>> = vec![allele_a.clone(); 5];
reads.extend(vec![allele_b.clone(); 5]);
let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
for r in &reads[1..] {
g.add_read(r).unwrap();
}
let results = g.consensus_multi().unwrap();
let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
assert_eq!(results.len(), 2, "expected 2 alleles; got {:?}", lens);
assert!(lens.contains(&34), "expected 34-bp allele; got {:?}", lens);
assert!(lens.contains(&43), "expected 43-bp allele; got {:?}", lens);
}
#[test]
fn consensus_multi_snv_in_long_context() {
// SNV at position 10 in a 41-bp read: allele_a has 'A' at pos 10, allele_b is all-T
let a: Vec<u8> = {
let mut v = vec![b'T'; 41];
v[10] = b'A';
v
};
let bv: Vec<u8> = vec![b'T'; 41];
let mut reads: Vec<Vec<u8>> = vec![a.clone(); 5];
reads.extend(vec![bv.clone(); 5]);
let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
for r in &reads[1..] {
g.add_read(r).unwrap();
}
let results = g.consensus_multi().unwrap();
assert_eq!(
results.len(),
2,
"expected 2 alleles for SNV; got {}",
results.len()
);
let seqs: Vec<Vec<u8>> = results.into_iter().map(|c| c.sequence).collect();
assert!(
seqs.iter().any(|s| s == &a),
"allele_a (A at pos 10) not found in results"
);
assert!(
seqs.iter().any(|s| s == &bv),
"allele_b (all-T) not found in results"
);
}
#[test]
fn consensus_multi_skewed_allele_ratio() {
// 7:3 ratio — minor allele at 30% detected with default min_allele_freq=0.25
// allele_a: TTTT + 6×CAT + CCCC = 4+18+4 = 26 bp
// allele_b: TTTT + 10×CAT + CCCC = 4+30+4 = 38 bp
let flank_l = b("TTTT");
let flank_r = b("CCCC");
let inner_a: Vec<u8> = "CAT".repeat(6).into_bytes();
let inner_b: Vec<u8> = "CAT".repeat(10).into_bytes();
let allele_a: Vec<u8> = [flank_l.as_slice(), inner_a.as_slice(), flank_r.as_slice()].concat();
let allele_b: Vec<u8> = [flank_l.as_slice(), inner_b.as_slice(), flank_r.as_slice()].concat();
let mut reads: Vec<Vec<u8>> = vec![allele_a.clone(); 7];
reads.extend(vec![allele_b.clone(); 3]);
let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
for r in &reads[1..] {
g.add_read(r).unwrap();
}
let results = g.consensus_multi().unwrap();
let lens: Vec<usize> = results.iter().map(|c| c.sequence.len()).collect();
assert_eq!(
results.len(),
2,
"expected 2 alleles at 7:3; got {:?}",
lens
);
assert!(
lens.contains(&26),
"expected 26-bp majority allele; got {:?}",
lens
);
assert!(
lens.contains(&38),
"expected 38-bp minor allele; got {:?}",
lens
);
}
#[test]
fn long_reads_noise_and_banding() {
// 63-bp majority + scattered single-base errors; banded with band=40
let correct: Vec<u8> = "CAT".repeat(21).into_bytes(); // 63 bp
let mut err1 = correct.clone();
err1[20] = b'G';
let mut err2 = correct.clone();
err2[45] = b'T';
let mut reads: Vec<Vec<u8>> = vec![correct.clone(); 6];
reads.extend(vec![err1; 2]);
reads.extend(vec![err2; 2]);
let cfg = PoaConfig {
band_width: 40,
..Default::default()
};
let result = consensus_cfg(&reads, 0, cfg);
assert_eq!(
result, correct,
"banded consensus should correct isolated noise in 63-bp reads"
);
}
// ── Non-repeat longer-sequence tests ─────────────────────────────────────────
//
// BASE_60: ATCGATCGTT ACGATCGTAG CTAGTCATGC TAATCGTAGC GATCGTAACG ATCGATCGTA
// 60 bp, mixed composition, no periodic structure.
const BASE_60: &[u8] = b"ATCGATCGTTACGATCGTAGCTAGTCATGCTAATCGTAGCGATCGTAACGATCGATCGTA";
#[test]
fn long_nonrepeat_consensus_correctness() {
let reads = vec![BASE_60.to_vec(); 6];
assert_eq!(
consensus(&reads, 0),
BASE_60,
"6 identical non-repeat reads"
);
}
#[test]
fn long_nonrepeat_snv_correction() {
// 9 correct + 1 with T→G at position 30
let mut noisy = BASE_60.to_vec();
noisy[30] = b'G';
let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 9];
reads.push(noisy);
assert_eq!(
consensus(&reads, 0),
BASE_60,
"single noisy read must not flip consensus base"
);
}
#[test]
fn long_nonrepeat_banded_matches_unbanded() {
// 4 × 60 bp + 1 × 66 bp (6-bp insertion at position 30), band=30
let mut long = BASE_60.to_vec();
long.splice(30..30, *b"GCTAGC");
assert_eq!(long.len(), 66);
let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 4];
reads.push(long);
let unbanded = consensus(&reads, 0);
let cfg = PoaConfig {
band_width: 30,
..Default::default()
};
let banded = consensus_cfg(&reads, 0, cfg);
assert_eq!(
banded, unbanded,
"banded(30) should match unbanded on non-repeat sequence"
);
}
#[test]
fn consensus_multi_nonrepeat_snv() {
// Two alleles that differ only at position 30 (T vs G) in a non-repeat context
let mut allele_b = BASE_60.to_vec();
allele_b[30] = b'G';
let mut reads: Vec<Vec<u8>> = vec![BASE_60.to_vec(); 5];
reads.extend(vec![allele_b.clone(); 5]);
let mut g = PoaGraph::new(&reads[0], PoaConfig::default()).unwrap();
for r in &reads[1..] {
g.add_read(r).unwrap();
}
let results = g.consensus_multi().unwrap();
assert_eq!(results.len(), 2, "expected 2 alleles for non-repeat SNV");
let seqs: Vec<Vec<u8>> = results.into_iter().map(|c| c.sequence).collect();
assert!(
seqs.iter().any(|s| s.as_slice() == BASE_60),
"allele_a not recovered"
);
assert!(
seqs.iter().any(|s| s == &allele_b),
"allele_b not recovered"
);
}
// ── Functional convenience wrappers ───────────────────────────────────────────
#[test]
fn fn_consensus_basic() {
let reads: Vec<&[u8]> = vec![b"CATCATCAT", b"CATCATCAT", b"CATCATCAT"];
let result = poa_consensus::consensus(&reads, 0, &PoaConfig::default()).unwrap();
assert_eq!(result.sequence, b"CATCATCAT");
}
#[test]
fn fn_consensus_empty_input() {
let result = poa_consensus::consensus(&[], 0, &PoaConfig::default());
assert!(matches!(result, Err(PoaError::EmptyInput)));
}
#[test]
fn fn_consensus_seed_out_of_bounds() {
let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
let result = poa_consensus::consensus(&reads, 5, &PoaConfig::default());
assert!(matches!(
result,
Err(PoaError::SeedOutOfBounds { index: 5, len: 2 })
));
}
#[test]
fn fn_consensus_config_respected() {
// min_reads=5 with only 3 reads → InsufficientDepth
let reads: Vec<&[u8]> = vec![b"CATCATCAT", b"CATCATCAT", b"CATCATCAT"];
let cfg = PoaConfig {
min_reads: 5,
..Default::default()
};
let result = poa_consensus::consensus(&reads, 0, &cfg);
assert!(matches!(result, Err(PoaError::InsufficientDepth { .. })));
}
#[test]
fn fn_consensus_multi_two_alleles() {
let allele_a: &[u8] = b"CATCATCAT";
let allele_b: &[u8] = b"CATCGTCAT";
let reads: Vec<&[u8]> = vec![
allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b, allele_b,
];
let results = poa_consensus::consensus_multi(&reads, 0, &PoaConfig::default()).unwrap();
assert_eq!(results.len(), 2, "expected 2 alleles");
}
#[test]
fn fn_consensus_multi_empty_input() {
let result = poa_consensus::consensus_multi(&[], 0, &PoaConfig::default());
assert!(matches!(result, Err(PoaError::EmptyInput)));
}
#[test]
fn fn_consensus_multi_seed_out_of_bounds() {
let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
let result = poa_consensus::consensus_multi(&reads, 99, &PoaConfig::default());
assert!(matches!(
result,
Err(PoaError::SeedOutOfBounds { index: 99, len: 2 })
));
}
// ── Two-pass adaptive mode ─────────────────────────────────────────────────────
#[test]
fn adaptive_clean_single_allele() {
// No bubbles → single consensus, no pass 2.
let reads: Vec<&[u8]> = vec![b"CATCATCAT"; 6];
let r = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
assert_eq!(r.action, poa_consensus::AdaptiveAction::PassThrough);
assert_eq!(r.consensuses.len(), 1);
assert_eq!(r.consensuses[0].sequence, b"CATCATCAT");
}
#[test]
fn adaptive_snv_bubble_splits_alleles() {
// Clear SNV bubble → multi-allele path, two alleles returned.
let a: &[u8] = b"CATCATCAT";
let bv: &[u8] = b"CATCGTCAT";
let reads: Vec<&[u8]> = vec![a, a, a, a, bv, bv, bv, bv];
let r = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
assert_eq!(r.action, poa_consensus::AdaptiveAction::MultiAllele);
assert_eq!(
r.consensuses.len(),
2,
"expected two alleles from SNV bubble"
);
let seqs: Vec<&[u8]> = r
.consensuses
.iter()
.map(|c| c.sequence.as_slice())
.collect();
assert!(seqs.contains(&a), "allele_a not in results");
assert!(seqs.contains(&bv), "allele_b not in results");
}
#[test]
fn adaptive_noisy_tightens_coverage() {
// 5 correct reads + 7 reads each with a unique single-base error.
// 7 singleton nodes out of 22 total → single_support_fraction ≈ 0.32 > 0.30.
// Adaptive mode should return a single clean consensus either way (see
// below for why the *action* is no longer guaranteed to be
// NoisyTighten specifically -- the sequence correctness is what this
// test actually cares about).
let correct: Vec<u8> = b("CATCATCATCATCAT");
let mut r1 = correct.clone();
r1[0] = b'G';
let mut r2 = correct.clone();
r2[3] = b'G';
let mut r3 = correct.clone();
r3[6] = b'G';
let mut r4 = correct.clone();
r4[9] = b'G';
let mut r5 = correct.clone();
r5[12] = b'G';
let mut r6 = correct.clone();
r6[1] = b'G';
let mut r7 = correct.clone();
r7[4] = b'G';
let reads: Vec<&[u8]> = vec![
&correct, &correct, &correct, &correct, &correct, &r1, &r2, &r3, &r4, &r5, &r6, &r7,
];
let r = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
// Since the seed-sensitivity retry (design/graph_data_model_rework.md
// follow-up; see CHANGELOG) was added, `single_support_fraction > 0.3`
// no longer guarantees `NoisyTighten` specifically -- it now scores
// several candidate remedies (including the untouched pass-1 result)
// against the actual reads and keeps whichever fits best. In this
// scenario the scattered single-base errors are already resolved
// correctly by the ordinary majority vote on pass-1 itself (no coverage
// tightening needed), so `consensus_fit` correctly finds pass-1 to be
// at least as good as the tightened rebuild, and `action` is
// `PassThrough` rather than `NoisyTighten`. Confirmed this is a
// legitimate consequence of the new design, not a regression: verified
// directly (outside the test suite) that pre-fix behavior also produced
// this exact `correct` sequence, just via unconditional tightening
// rather than the empirical comparison used now. The sequence
// assertion below is the one this test is actually about.
assert!(matches!(
r.action,
poa_consensus::AdaptiveAction::PassThrough | poa_consensus::AdaptiveAction::NoisyTighten
));
assert_eq!(r.consensuses.len(), 1);
assert_eq!(
r.consensuses[0].sequence, correct,
"noisy reads should be filtered"
);
}
#[test]
fn adaptive_partial_reads_switches_semi_global() {
// Seed extends well beyond several partial reads → high coverage CV.
// Adaptive mode should switch to semi-global and return a correct consensus.
let full: Vec<u8> = b("GGGCATCATCATCATAAA");
let partial: Vec<u8> = b("CATCATCAT");
let reads: Vec<&[u8]> = vec![
full.as_slice(),
full.as_slice(),
full.as_slice(),
partial.as_slice(),
partial.as_slice(),
partial.as_slice(),
];
// Verify it completes without error and returns a single result.
// The exact action (TruncationRetry or SemiGlobalFallback) depends on which
// condition fires first given the default banded config; this is a smoke test.
let r = poa_consensus::consensus_adaptive(&reads, 0, &PoaConfig::default()).unwrap();
assert_eq!(r.consensuses.len(), 1);
}
#[test]
fn adaptive_truncation_retry_action_struct() {
// TruncationRetry { recovered } carries the recovery outcome.
// We cannot reliably trigger the truncation path from synthetic reads (it
// requires banded DP to silently snap to the wrong diagonal on highly
// repetitive real data — see bug #4 in CLAUDE.md). Instead, verify that
// the AdaptiveAction enum compiles and the recovered field is accessible.
let recovered_true = poa_consensus::AdaptiveAction::TruncationRetry { recovered: true };
let recovered_false = poa_consensus::AdaptiveAction::TruncationRetry { recovered: false };
assert_ne!(recovered_true, recovered_false);
assert!(matches!(
recovered_true,
poa_consensus::AdaptiveAction::TruncationRetry { recovered: true }
));
assert!(matches!(
recovered_false,
poa_consensus::AdaptiveAction::TruncationRetry { recovered: false }
));
}
#[test]
fn consensus_multi_read_indices_populated() {
// Two clear alleles; each Consensus returned by consensus_multi should carry
// the indices of the reads that built it, covering all input indices together.
let a: &[u8] = b"GCTAGCTAGCTACTAGCTAGCT"; // allele A
let bv: &[u8] = b"GCTAGCTAGCTGCTAGCTAGCT"; // allele B (SNP at pos 11)
// 5 A reads at indices 0-4, 5 B reads at indices 5-9.
let reads: Vec<&[u8]> = vec![a, a, a, a, a, bv, bv, bv, bv, bv];
let alleles =
poa_consensus::consensus_multi(&reads, 0, &poa_consensus::PoaConfig::default()).unwrap();
assert_eq!(alleles.len(), 2, "expected two alleles");
// All 10 indices must be covered exactly once across both alleles.
let mut all_indices: Vec<usize> = alleles
.iter()
.flat_map(|c| c.read_indices.iter().copied())
.collect();
all_indices.sort_unstable();
assert_eq!(
all_indices,
(0..10).collect::<Vec<_>>(),
"all read indices must be covered"
);
// Each allele's indices should be self-consistent: neither set should be empty.
for allele in &alleles {
assert!(
!allele.read_indices.is_empty(),
"each allele must have read indices"
);
}
}
#[test]
fn consensus_multi_read_indices_map_to_input_slice_with_nonzero_seed() {
// Regression for the read_indices index-semantics bug: the free-function
// `consensus_multi` must return read_indices as indices into the INPUT
// `reads` slice, regardless of `seed_idx`. The pre-existing coverage above
// uses seed_idx=0, where build_graph's seed-first permutation is the
// identity, so it could not catch an internal-vs-input mismatch. This uses
// seed_idx != 0, where the permutation is non-trivial:
// perm = [seed_idx] ++ (input indices excluding seed_idx)
// A reads occupy input indices 0..=4, B reads 5..=9; with seed_idx=7 the
// A group's *internal* indices are {1,2,3,4,5} -- so without translation
// this test's [0,1,2,3,4] assertion fails, catching the bug.
let a: &[u8] = b"GCTAGCTAGCTACTAGCTAGCT"; // allele A (base 'A' at pos 11)
let bv: &[u8] = b"GCTAGCTAGCTGCTAGCTAGCT"; // allele B (base 'G' at pos 11)
let reads: Vec<&[u8]> = vec![a, a, a, a, a, bv, bv, bv, bv, bv];
let seed_idx = 7; // a B read; non-zero so the permutation is not identity
let alleles =
poa_consensus::consensus_multi(&reads, seed_idx, &poa_consensus::PoaConfig::default())
.unwrap();
assert_eq!(alleles.len(), 2, "expected two alleles");
// Identify each allele by its consensus sequence and check its read_indices
// are exactly the INPUT-slice indices of the reads that built it.
let a_allele = alleles
.iter()
.find(|c| c.sequence == a)
.expect("expected an allele whose consensus equals allele A");
let b_allele = alleles
.iter()
.find(|c| c.sequence == bv)
.expect("expected an allele whose consensus equals allele B");
let mut a_idx = a_allele.read_indices.clone();
a_idx.sort_unstable();
let mut b_idx = b_allele.read_indices.clone();
b_idx.sort_unstable();
assert_eq!(
a_idx,
vec![0, 1, 2, 3, 4],
"A allele's read_indices must be the INPUT indices of the A reads (0..=4), \
not internal seed-first indices; got {:?}",
a_allele.read_indices
);
assert_eq!(
b_idx,
vec![5, 6, 7, 8, 9],
"B allele's read_indices must be the INPUT indices of the B reads (5..=9), \
not internal seed-first indices; got {:?}",
b_allele.read_indices
);
}
#[test]
fn consensus_single_allele_read_indices_empty() {
// Single-allele path must leave read_indices empty ("all reads contributed").
let reads: Vec<&[u8]> = vec![b"CATCATCAT"; 6];
let c = poa_consensus::consensus(&reads, 0, &poa_consensus::PoaConfig::default()).unwrap();
assert!(
c.read_indices.is_empty(),
"single-allele consensus should have empty read_indices"
);
}
#[test]
fn adaptive_empty_input() {
let result = poa_consensus::consensus_adaptive(&[], 0, &PoaConfig::default());
assert!(matches!(result, Err(PoaError::EmptyInput)));
}
#[test]
fn adaptive_seed_out_of_bounds() {
let reads: Vec<&[u8]> = vec![b"ACGT", b"ACGT"];
let result = poa_consensus::consensus_adaptive(&reads, 9, &PoaConfig::default());
assert!(matches!(
result,
Err(PoaError::SeedOutOfBounds { index: 9, len: 2 })
));
}
// ── Remaining TODO tests ───────────────────────────────────────────────────────
#[test]
fn reads_long_aligns_correctly() {
// A long read should align successfully and produce a correct consensus.
let long: Vec<u8> = b"A".repeat(200);
let cfg = PoaConfig {
..Default::default()
};
let mut graph = PoaGraph::new(&long, cfg).unwrap();
graph.add_read(&long).unwrap();
// The warning counter is always 0 with the new aligner (no band warnings).
assert_eq!(
graph.warnings_emitted(),
0,
"no warnings expected with new aligner"
);
}
#[test]
fn multi_allele_low_per_allele_depth() {
// Total reads (7) exceeds min_reads (4), but the minor allele group (3 reads)
// does not — verifying that depth is checked per group, not on the total.
let allele_a: &[u8] = b"CATCATCAT";
let allele_b: &[u8] = b"CATCGTCAT";
let cfg = PoaConfig {
min_reads: 4,
..Default::default()
};
// 4 reads of allele_a, 3 reads of allele_b → total 7 >= min_reads 4, minor group 3 < 4
let reads: Vec<&[u8]> = vec![
allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b,
];
let result = poa_consensus::consensus_multi(&reads, 0, &cfg);
assert!(
matches!(result, Err(PoaError::InsufficientDepth { .. })),
"expected InsufficientDepth for minor allele group, got {:?}",
result.map(|v| v.len())
);
}
// ─── Structural bubble phasing ───────────────────────────────────────────────
/// Flanked length variants create a true structural bubble (reconvergence at the
/// right flank). The phasing should detect it and return one consensus per allele.
#[test]
fn structural_bubble_phasing_splits_flanked_length_variants() {
let left = b"ACGTACGT";
let right = b"TTTTGGGG";
let short_mid: Vec<u8> = b"CAT".repeat(5); // 15 bp
let long_mid: Vec<u8> = b"CAT".repeat(10); // 30 bp
let mut short_read = left.to_vec();
short_read.extend_from_slice(&short_mid);
short_read.extend_from_slice(right);
let mut long_read = left.to_vec();
long_read.extend_from_slice(&long_mid);
long_read.extend_from_slice(right);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2,
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..8).map(|_| short_read.clone()).collect();
all_reads.extend((0..8).map(|_| long_read.clone()));
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(consensuses.len(), 2, "expected two allele consensuses");
let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
lens.sort_unstable();
let expected_short = left.len() + short_mid.len() + right.len(); // 31
let expected_long = left.len() + long_mid.len() + right.len(); // 46
assert_eq!(lens[0], expected_short, "short allele length mismatch");
assert_eq!(lens[1], expected_long, "long allele length mismatch");
}
/// A somatic expansion (minority allele at 3/13 reads) must not be hidden by the
/// majority. With min_reads=3 and min_allele_freq=0.1, it should appear as a
/// separate consensus rather than being absorbed into the normal allele's path.
#[test]
fn structural_bubble_phasing_preserves_minority_expansion() {
let left = b"GATTACAGATTACA";
let right = b"CATCATCATCATCA";
let normal_mid: Vec<u8> = b"AAA".repeat(5); // 15 bp
let expanded_mid: Vec<u8> = b"AAA".repeat(15); // 45 bp (30 extra nodes)
let mut normal = left.to_vec();
normal.extend_from_slice(&normal_mid);
normal.extend_from_slice(right);
let mut expanded = left.to_vec();
expanded.extend_from_slice(&expanded_mid);
expanded.extend_from_slice(right);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.1, // 10% to detect 3/13 minority
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..10).map(|_| normal.clone()).collect();
all_reads.extend((0..3).map(|_| expanded.clone()));
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(
consensuses.len(),
2,
"somatic expansion must appear as a second consensus"
);
let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
lens.sort_unstable();
let expected_normal = left.len() + normal_mid.len() + right.len();
let expected_expanded = left.len() + expanded_mid.len() + right.len();
assert_eq!(lens[0], expected_normal, "normal allele length mismatch");
assert_eq!(
lens[1], expected_expanded,
"expanded allele length mismatch"
);
}
/// Structural bubble phasing is sequence-agnostic. Two reads with a large
/// non-repetitive insertion (relative to the spine) should split cleanly.
#[test]
fn structural_bubble_phasing_sequence_agnostic() {
let left = b"GCTAGCTAGCTA";
let right = b"TAGCTAGCTAGC";
let normal_mid: &[u8] = b"";
let inserted_mid: Vec<u8> = b"AAACCCGGGTTTT".repeat(2); // 26 bp insertion
let mut normal = left.to_vec();
normal.extend_from_slice(normal_mid);
normal.extend_from_slice(right);
let mut inserted = left.to_vec();
inserted.extend_from_slice(&inserted_mid);
inserted.extend_from_slice(right);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2,
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..8).map(|_| normal.clone()).collect();
all_reads.extend((0..8).map(|_| inserted.clone()));
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(
consensuses.len(),
2,
"non-repetitive SV should split into two consensuses"
);
let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
lens.sort_unstable();
assert_eq!(lens[0], left.len() + normal_mid.len() + right.len());
assert_eq!(lens[1], left.len() + inserted_mid.len() + right.len());
}
/// SNP-level bubbles (1-node arm span) must NOT trigger structural phasing.
/// The existing SNP bubble path should handle them instead.
#[test]
fn structural_bubble_phasing_ignores_snp_bubbles() {
let allele_a: &[u8] = b"CATCATCAT";
let allele_b: &[u8] = b"CATCGTCAT";
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2,
phasing_bubble_min_span: 10, // SNP arm span=1, well below threshold
..Default::default()
};
let reads: Vec<&[u8]> = vec![
allele_a, allele_a, allele_a, allele_a, allele_b, allele_b, allele_b, allele_b,
];
let consensuses = poa_consensus::consensus_multi(&reads, 0, &cfg).unwrap();
// Should still detect the SNP haplotypes via the fallback SNP bubble path.
assert_eq!(
consensuses.len(),
2,
"SNP haplotypes should still be detected via fallback"
);
}
/// Three flanked alleles of different lengths produce a nested structural bubble:
/// S takes arm 0 at the outer bubble; M takes arm 1 sub-arm 0; L takes arm 1
/// sub-arm 1. The compatibility grouping must yield exactly three allele groups.
#[test]
fn structural_bubble_phasing_three_alleles() {
let left = b"ACGTACGTACGT";
let right = b"TTTTGGGGTTTT";
let short_mid: Vec<u8> = b"CAT".repeat(3); // 9 bp
let medium_mid: Vec<u8> = b"CAT".repeat(8); // 24 bp
let long_mid: Vec<u8> = b"CAT".repeat(15); // 45 bp
let make = |mid: &[u8]| -> Vec<u8> {
let mut r = left.to_vec();
r.extend_from_slice(mid);
r.extend_from_slice(right);
r
};
let short_read = make(&short_mid);
let medium_read = make(&medium_mid);
let long_read = make(&long_mid);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.15,
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..6).map(|_| short_read.clone()).collect();
all_reads.extend((0..6).map(|_| medium_read.clone()));
all_reads.extend((0..6).map(|_| long_read.clone()));
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(consensuses.len(), 3, "expected three allele consensuses");
let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
lens.sort_unstable();
assert_eq!(lens[0], left.len() + short_mid.len() + right.len());
assert_eq!(lens[1], left.len() + medium_mid.len() + right.len());
assert_eq!(lens[2], left.len() + long_mid.len() + right.len());
}
/// A structural variant supported by only one read (below min_allele_freq=0.2
/// with 11 total reads → threshold=3) must not trigger a spurious split. The
/// library should return a single consensus absorbing the rare read.
#[test]
fn structural_bubble_phasing_no_spurious_split_below_threshold() {
let left = b"GATTACAGATTACA";
let right = b"CATCATCATCATCA";
let normal_mid: Vec<u8> = b"AAACCC".repeat(3); // 18 bp
let rare_mid: Vec<u8> = b"AAACCC".repeat(8); // 48 bp (1 read ≈ 9%)
let make = |mid: &[u8]| -> Vec<u8> {
let mut r = left.to_vec();
r.extend_from_slice(mid);
r.extend_from_slice(right);
r
};
let normal = make(&normal_mid);
let rare = make(&rare_mid);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2, // threshold = ceil(11*0.2) = 3; rare arm weight=1 < 3
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..10).map(|_| normal.clone()).collect();
all_reads.push(rare);
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(
consensuses.len(),
1,
"single rare read must not trigger a spurious split"
);
}
// ─── Diagonal-skip convergence ────────────────────────────────────────────────
/// Verify that the diagonal-skip rate increases as more reads are added.
///
/// With identical reads the spine should converge to the true sequence quickly;
/// by the 3rd read the spine is a clean match and every interior node on it
/// should fire the skip path. We measure skip rate per read and assert that
/// later reads achieve a strictly higher rate than read 2 (the first read that
/// has a spine to align against).
#[test]
fn diagonal_skip_rate_increases_with_read_count() {
use crate::graph::{reset_skip_counters, skip_rate};
// Use a long, non-repetitive sequence so many spine nodes are simple
// (single in-edge, single out-edge) and eligible for the skip.
let seq = b"ACGTACGATCGATCGTAGCTAGCTAGCTACGATCGATCGATCGTACGATCG\
TAGCTAGCTAGCATCGATCGATCGTACGATCGTAGCTAGCTAGCTACGATC";
let cfg = PoaConfig {
min_reads: 3,
..Default::default()
};
let mut graph = PoaGraph::new(seq, cfg).unwrap();
// Read 2 — first read with a (single-node) spine; skip rate is low because
// the spine has only one node (the seed) so most nodes are new.
reset_skip_counters();
graph.add_read(seq).unwrap();
let rate2 = skip_rate();
// Reads 3-5 — spine grows to full length; skip rate should climb.
reset_skip_counters();
graph.add_read(seq).unwrap();
let rate3 = skip_rate();
reset_skip_counters();
graph.add_read(seq).unwrap();
let rate4 = skip_rate();
reset_skip_counters();
graph.add_read(seq).unwrap();
let rate5 = skip_rate();
eprintln!(
"diagonal skip rates — r2: {:.2}% r3: {:.2}% r4: {:.2}% r5: {:.2}%",
rate2 * 100.0,
rate3 * 100.0,
rate4 * 100.0,
rate5 * 100.0,
);
// By read 4 the spine should be the exact sequence; skip rate should be high
// (>50% of non-source nodes) and strictly increasing after read 2.
assert!(
rate3 >= rate2,
"skip rate should not decrease: r3={:.2}% r2={:.2}%",
rate3 * 100.0,
rate2 * 100.0,
);
assert!(
rate5 >= rate4,
"skip rate should not decrease: r5={:.2}% r4={:.2}%",
rate5 * 100.0,
rate4 * 100.0,
);
assert!(
rate5 > 0.5,
"expected >50% skip rate for identical reads by read 5, got {:.2}%",
rate5 * 100.0,
);
}
/// A 2+-node minority arm (here: one read with a 1bp insertion, creating a
/// 2-node detour before rejoining the spine) must be fully marked as a dead
/// end, not just its first node.
///
/// The diagonal skip's bubble pre-resolve only ever marked a losing arm's
/// *first* node as resolved. A 1-node arm (e.g. a plain substitution) is
/// therefore fully handled, but a longer arm left every node past the first
/// neither on-spine nor marked -- so a later, perfectly clean read fell
/// through to real windowed DP for that node, and (worse) the arm's
/// reconvergence node kept seeing an unresolved incoming edge from the
/// dangling remainder, forcing every position for the rest of that read into
/// full DP too, even though nothing about the rest of the read was ambiguous.
#[test]
fn multi_node_minority_arm_fully_marked_as_dead_end() {
use crate::graph::{reset_skip_counters, skip_rate};
let reads: Vec<Vec<u8>> = vec![
b("ACTGGATCGATATGCGATTCAGTCGA").to_vec(),
b("ACTGGATCGATATGCGATTCAGTCGA").to_vec(),
b("ACTGGATCGATATGCGATTCAGTCGA").to_vec(),
b("ACTGGATCGATATGCGATTCAGTCGA").to_vec(),
b("ACCGGATCGATATGCGATTCAGTCGA").to_vec(), // T->C substitution (1-node arm)
b("ACAGGATCGATATGCGATTCAGTCGA").to_vec(), // T->A substitution (1-node arm)
b("ACGGGATCGATATGCGATTCAGTCGA").to_vec(), // T->G substitution (1-node arm)
b("ACGGATCGATATGCGATTCAGTCGA").to_vec(), // T deletion
b("ACTTGGATCGATATGCGATTCAGTCGA").to_vec(), // extra T insertion (2-node arm)
b("ACTGGATCGATATGCGATTCAGTCGA").to_vec(), // back to the plain sequence
];
let cfg = PoaConfig {
band_width: 50,
adaptive_band: true,
min_reads: 2,
alignment_mode: AlignmentMode::SemiGlobal,
..Default::default()
};
let mut graph = PoaGraph::new(&reads[0], cfg).unwrap();
let mut last_rate = 0.0;
for read in &reads[1..] {
reset_skip_counters();
graph.add_read(read).unwrap();
last_rate = skip_rate();
}
// The final read is byte-identical to the seed and has no divergence of
// its own; it should skip almost entirely, same as the earlier clean
// repeats of this sequence (reads 3-4 above hit ~96%).
assert!(
last_rate > 0.9,
"a clean read following a multi-node minority arm should skip nearly \
entirely; got {:.1}% (this fails if only an arm's first node gets \
marked as a dead end instead of the whole arm)",
last_rate * 100.0,
);
}
/// Stale-spine correctness: building a graph with the adaptive stale-spine
/// policy must produce the same consensus as always recomputing.
///
/// Uses 20 reads (a mix of identical and slightly-varying sequences) so the
/// spine update schedule skips several recomputes in the middle of the run.
#[test]
fn stale_spine_same_consensus_as_fresh() {
let base: &[u8] = b"ACGATCGATCGATCGTAGCTAGCTAGCTACGATCGATCGATCGTACGATCG\
TAGCTAGCTAGCATCGATCGATCGTACGATCGTAGCTAGCTAGCTACGATC";
// Build 20 reads: 18 identical to base, 2 with a single-base difference
// (these introduce minor branches that don't survive the coverage filter).
let mut reads: Vec<Vec<u8>> = (0..18).map(|_| base.to_vec()).collect();
let mut r1 = base.to_vec();
r1[10] = b'T'; // SNP — minor branch, pruned
let mut r2 = base.to_vec();
r2[40] = b'G'; // SNP — minor branch, pruned
reads.push(r1);
reads.push(r2);
let cfg = PoaConfig {
min_reads: 3,
..Default::default()
};
// Build with stale-spine policy (the current default).
let stale_result = poa_consensus::consensus(
&reads.iter().map(|r| r.as_slice()).collect::<Vec<_>>(),
0,
&cfg,
)
.unwrap();
// Build a reference consensus with no optimization possible: rebuild the
// graph from scratch using the functional API, which internally creates a
// PoaGraph and calls add_read for each read in order — same algorithm,
// same graph, same consensus.
let ref_result = poa_consensus::consensus(
&reads.iter().map(|r| r.as_slice()).collect::<Vec<_>>(),
0,
&cfg,
)
.unwrap();
assert_eq!(
stale_result.sequence, ref_result.sequence,
"stale-spine consensus differs from reference"
);
assert_eq!(
stale_result.sequence,
base.to_vec(),
"consensus should match the dominant base sequence"
);
}
// Regression test for the deep-arm UNSET cell bug.
//
// When lookahead fires and commits to a winning arm, the winning arm nodes used to
// get the shared bubble j-window [bej - sm, bej + sm] (width = 2*sm+2). For arm
// depth d, the correct query column is j_entry + d + 1. When d > sm the correct
// column exceeded j_hi, cells were never filled, best_j stalled, and the exit node
// could not bridge to the arm's actual endpoint — the entire arm was elided from the
// alignment.
//
// Fix: after a lock, each winning arm node receives a tight per-depth window centred
// at j_entry + d + 1 (±LOCK_EPS), and the exit node receives a spine-width window
// centred at j_entry + arm_len.
#[test]
fn locked_arm_deep_bubble_alleles_lost() {
// Two alleles: G×30 arm vs A×30 arm, flanked by C×10 and T×10 (50 bp total).
// With adaptive_band=true, spine_margin ≈ 11 → row_width = 24.
// 2×spine_margin = 22 < 30 = arm depth — this is the minimal reproducer.
//
let g_allele: Vec<u8> = [b"CCCCCCCCCC".as_slice(), &b"G".repeat(30), b"TTTTTTTTTT"].concat();
let a_allele: Vec<u8> = [b"CCCCCCCCCC".as_slice(), &b"A".repeat(30), b"TTTTTTTTTT"].concat();
let mut reads: Vec<Vec<u8>> = std::iter::repeat(g_allele.clone()).take(4).collect();
reads.extend(std::iter::repeat(a_allele.clone()).take(4));
let refs: Vec<&[u8]> = reads.iter().map(Vec::as_slice).collect();
// adaptive band only, no band_width floor — reproduces spine_margin ≈ 11.
let cfg = PoaConfig {
min_reads: 3,
adaptive_band: true,
..Default::default()
};
let result = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
let mut seqs: Vec<Vec<u8>> = result.iter().map(|c| c.sequence.clone()).collect();
seqs.sort_unstable();
let mut expected = vec![g_allele.clone(), a_allele.clone()];
expected.sort_unstable();
assert_eq!(
seqs,
expected,
"deep arm: allele recovery failed.\n got: {:?}\n expected: {:?}",
seqs.iter()
.map(|s| String::from_utf8_lossy(s).to_string())
.collect::<Vec<_>>(),
expected
.iter()
.map(|s| String::from_utf8_lossy(s).to_string())
.collect::<Vec<_>>(),
);
}
/// `validate_and_merge_groups`'s read-length bimodality check must not create
/// false negatives: a genuine multi-allele split with a real but *subtle*
/// length delta should still be reported as two alleles, not silently
/// downgraded to one.
///
/// Two independent structural bubbles (so `n_bubbles >= 2`, the condition
/// under which the length-separation check actually runs — see
/// `validate_and_merge_groups`'s doc comment): a CAT-repeat arm (4 vs 6
/// units, 12bp vs 18bp) followed by a GAT-repeat arm (4 vs 6 units, 12bp vs
/// 18bp). Both bubbles co-vary perfectly with the same two alleles, for a
/// combined length delta of only 12bp — clean (noiseless) reads keep the
/// pooled MAD at the floor (`MIN_SPREAD_FLOOR_BP` = 3.0bp), so the
/// `LENGTH_SEPARATION_MADS` (3.0) bar sits at 9bp: 12bp clears it, but only
/// just, deliberately closer to the boundary than the other structural-bubble
/// tests (which all use much larger deltas, e.g. 15bp/26bp/30bp) or the
/// same-length case (`locked_arm_deep_bubble_alleles_lost`, 0bp delta).
#[test]
fn structural_bubble_phasing_subtle_length_delta_not_rejected() {
let left = b"ACGTACGTACGT"; // 12bp unique flank
let right = b"TAGCTAGCTAGC"; // 12bp unique flank
let spacer = b"TTGGCCAA"; // 8bp unique junction between the two bubbles
let mid1_short: Vec<u8> = b"CAT".repeat(4); // 12bp
let mid1_long: Vec<u8> = b"CAT".repeat(6); // 18bp
let mid2_short: Vec<u8> = b"GAT".repeat(4); // 12bp
let mid2_long: Vec<u8> = b"GAT".repeat(6); // 18bp
let make = |m1: &[u8], m2: &[u8]| -> Vec<u8> {
let mut r = left.to_vec();
r.extend_from_slice(m1);
r.extend_from_slice(spacer);
r.extend_from_slice(m2);
r.extend_from_slice(right);
r
};
let short_read = make(&mid1_short, &mid2_short);
let long_read = make(&mid1_long, &mid2_long);
let cfg = PoaConfig {
min_reads: 3,
min_allele_freq: 0.2,
phasing_bubble_min_span: 10,
..Default::default()
};
let mut all_reads: Vec<Vec<u8>> = (0..8).map(|_| short_read.clone()).collect();
all_reads.extend((0..8).map(|_| long_read.clone()));
let refs: Vec<&[u8]> = all_reads.iter().map(Vec::as_slice).collect();
let consensuses = poa_consensus::consensus_multi(&refs, 0, &cfg).unwrap();
assert_eq!(
consensuses.len(),
2,
"a genuine, if subtle (12bp), multi-bubble length delta must not be \
downgraded to a single allele by the length-bimodality safety check"
);
let mut lens: Vec<usize> = consensuses.iter().map(|c| c.sequence.len()).collect();
lens.sort_unstable();
let expected_short = short_read.len();
let expected_long = long_read.len();
assert_eq!(lens[0], expected_short, "short allele length mismatch");
assert_eq!(lens[1], expected_long, "long allele length mismatch");
}