use std::{
collections::VecDeque,
fmt::{self, Display},
ops::{Index, IndexMut},
time::{Duration, Instant},
};
use anyhow::anyhow;
use bstr::ByteSlice;
use compact_genome::interface::alphabet::AlphabetCharacter;
use generic_a_star::cost::AStarCost;
use lib_tsalign::{
a_star_aligner::{
alignment_geometry::AlignmentRange,
alignment_result::alignment::Alignment,
configurable_a_star_align::alphabets::dna_alphabet_or_n::{
DnaAlphabetOrN, DnaCharacterOrN,
},
template_switch_distance::{
AlignmentType, EqualCostRange, TemplateSwitchAncestor, TemplateSwitchDescendant,
TemplateSwitchDirection,
},
},
config::TemplateSwitchConfig,
costs::U64Cost,
};
use serde::{Deserialize, Serialize};
use tracing::instrument;
use crate::common::{
ImmutableSequence,
aligner::result::{
AlignmentFailure, AlignmentWithCost, TwitcherAlignmentResult,
TwitcherAlignmentWithStatistics,
},
};
#[derive(Deserialize, Serialize, Default)]
pub struct FourPointAligner {
costs: TemplateSwitchConfig<DnaAlphabetOrN, U64Cost>,
no_ts: bool,
}
struct AlignerInstance<'a> {
start_time: Instant,
reference: ImmutableSequence,
query: ImmutableSequence,
ranges: AlignmentRange,
costs: &'a TemplateSwitchConfig<DnaAlphabetOrN, U64Cost>,
no_ts: bool,
pre_ts: Mat<Cell>,
ts_qrr: Mat<TsCell>,
ts_rqr: Mat<TsCell>,
ts_qrf: Mat<TsCell>,
ts_rqf: Mat<TsCell>,
post_ts: Mat<Cell>,
}
#[allow(dead_code)]
#[derive(Debug, Clone, Deserialize, Serialize)]
pub struct FpaAlignmentStatistics {
pub duration: Duration,
pub ranges: AlignmentRange,
}
impl FourPointAligner {
pub const fn new(costs: TemplateSwitchConfig<DnaAlphabetOrN, U64Cost>, no_ts: bool) -> Self {
Self { costs, no_ts }
}
#[instrument(skip_all)]
pub fn align(
&self,
reference: ImmutableSequence,
query: ImmutableSequence,
ranges: AlignmentRange,
) -> TwitcherAlignmentResult {
let instance = AlignerInstance::new(reference, query, ranges, &self.costs, self.no_ts)?;
let res = instance.run()?;
Ok(res)
}
}
impl<'a> AlignerInstance<'a> {
fn new(
reference: ImmutableSequence,
query: ImmutableSequence,
ranges: AlignmentRange,
costs: &'a TemplateSwitchConfig<DnaAlphabetOrN, U64Cost>,
no_ts: bool,
) -> Result<Self, AlignmentFailure> {
let start_time = Instant::now();
// Forward alignment within the cluster range
let pre_ts: Mat<Cell> = Mat::try_new(
ranges.reference_range().len() + 1,
ranges.query_range().len() + 1,
)?;
// The TS matrices and post_ts are allocated lazily in `run()`, once `pre_ts` has been
// filled: a direction is only worth filling if it could possibly beat the no-TS
// alignment (see `run`). Start them empty.
Ok(Self {
start_time,
reference,
query,
ranges,
costs,
no_ts,
pre_ts,
ts_qrr: Mat::empty(),
ts_rqr: Mat::empty(),
ts_qrf: Mat::empty(),
ts_rqf: Mat::empty(),
post_ts: Mat::empty(),
})
}
fn run(mut self) -> Result<TwitcherAlignmentWithStatistics, AlignmentFailure> {
self.fill_pre_ts()?;
if !self.no_ts {
// Cost of the best no-TS alignment (pre_ts bottom-right corner). Any TS via
// direction `d` pays at least its base cost (everything else along the path is
// non-negative), so it can only beat the no-TS alignment when `base_cost_d` is
// strictly cheaper than this corner. Directions disabled via a MAX base cost are
// covered for free. Skipping a direction avoids both its allocation and its fill.
let pre_ts_corner = self
.pre_ts
.fin()
.map_or(U64Cost::from_primitive(u64::MAX), |c| c.score);
let mut any_ts = false;
// Fill each secondary matrix whose base cost could still beat the no-TS alignment.
// TODO since we only consider matches, surely the secondary fills can be done in
// linear space at least... :D
for kind in ALL_TS_KINDS {
let base = self
.costs
.base_cost
.get(kind.descendant, kind.ancestor, kind.direction);
if base < pre_ts_corner {
let (rows, cols) = self.ts_dims(kind);
let mut mat = Mat::try_new(rows, cols)?;
self.fill_ts(kind, &mut mat)?;
*self.ts_mat_mut(kind.which) = mat;
any_ts = true;
}
}
// post_ts (TS followed by more forward alignment) is reachable only via a TS jump,
// so it is pointless unless at least one direction is active.
if any_ts {
self.post_ts = Mat::try_new(
self.ranges.reference_range().len() + 1,
self.ranges.query_range().len() + 1,
)?;
self.fill_post_ts()?;
}
}
let ts_alignment = self.retrace_path()?;
let statistics = FpaAlignmentStatistics {
duration: self.start_time.elapsed(),
ranges: self.ranges,
};
Ok(TwitcherAlignmentWithStatistics {
alignment: ts_alignment,
stats: statistics.into(),
})
}
fn fill_pre_ts(&mut self) -> anyhow::Result<()> {
for i in 0..self.pre_ts.m() {
let rc = (i > 0)
.then(|| get_char(&self.reference, self.ranges.reference_offset() + i - 1))
.transpose()?;
for j in 0..self.pre_ts.n() {
let qc = (j > 0)
.then(|| get_char(&self.query, self.ranges.query_offset() + j - 1))
.transpose()?;
let matches = rc == qc;
if let (Some(rc), Some(qc)) = (rc, qc) {
// Mis-/Match
let cost = self
.costs
.primary_edit_costs
.match_or_substitution_cost(rc, qc);
let score = self.pre_ts[(i - 1, j - 1)].score + cost;
if self.pre_ts[(i, j)].ty.is_none() || score < self.pre_ts[(i, j)].score {
self.pre_ts[(i, j)].score = self.pre_ts[(i - 1, j - 1)].score + cost;
self.pre_ts[(i, j)].ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
}
}
if let Some(rc) = rc {
let cost = self.costs.primary_edit_costs.gap_extend_cost(rc);
let score = self.pre_ts[(i - 1, j)].score + cost;
if self.pre_ts[(i, j)].ty.is_none() || score < self.pre_ts[(i, j)].score {
self.pre_ts[(i, j)].score = score;
self.pre_ts[(i, j)].ty = Some(Align::Deletion);
}
}
if let Some(qc) = qc {
let cost = self.costs.primary_edit_costs.gap_extend_cost(qc);
let score = self.pre_ts[(i, j - 1)].score + cost;
if self.pre_ts[(i, j)].ty.is_none() || score < self.pre_ts[(i, j)].score {
self.pre_ts[(i, j)].score = score;
self.pre_ts[(i, j)].ty = Some(Align::Insertion);
}
}
}
}
Ok(())
}
/// Logical extents: full ancestor sequence length, descendant cluster-range length.
fn ts_extents(&self, kind: TsKind) -> (usize, usize) {
if kind.query_descendant() {
(self.reference.len(), self.ranges.query_range().len())
} else {
(self.query.len(), self.ranges.reference_range().len())
}
}
/// Physical `(rows, cols)` to allocate for `kind`.
fn ts_dims(&self, kind: TsKind) -> (usize, usize) {
let (a_len, d_len) = self.ts_extents(kind);
kind.phys(a_len + 1, d_len + 1)
}
fn ts_mat(&self, which: WhichMatrix) -> &Mat<TsCell> {
match which {
WhichMatrix::TsQrr => &self.ts_qrr,
WhichMatrix::TsRqr => &self.ts_rqr,
WhichMatrix::TsQrf => &self.ts_qrf,
WhichMatrix::TsRqf => &self.ts_rqf,
WhichMatrix::PreTs | WhichMatrix::PostTs => {
unreachable!("not a secondary matrix")
}
}
}
fn ts_mat_mut(&mut self, which: WhichMatrix) -> &mut Mat<TsCell> {
match which {
WhichMatrix::TsQrr => &mut self.ts_qrr,
WhichMatrix::TsRqr => &mut self.ts_rqr,
WhichMatrix::TsQrf => &mut self.ts_qrf,
WhichMatrix::TsRqf => &mut self.ts_rqf,
WhichMatrix::PreTs | WhichMatrix::PostTs => {
unreachable!("not a secondary matrix")
}
}
}
/// Ancestor character at logical index `a` (1-based, over the full ancestor sequence),
/// reverse-complemented for a reverse template switch.
fn ts_ancestor_char(&self, kind: TsKind, a: usize) -> anyhow::Result<DnaCharacterOrN> {
let seq = match kind.ancestor {
TemplateSwitchAncestor::Reference => &self.reference,
TemplateSwitchAncestor::Query => &self.query,
};
let c = get_char(seq, a - 1)?;
Ok(match kind.direction {
TemplateSwitchDirection::Reverse => c.complement(),
TemplateSwitchDirection::Forward => c,
})
}
/// Descendant character at logical index `d` (1-based, offset into the cluster range).
fn ts_descendant_char(&self, kind: TsKind, d: usize) -> anyhow::Result<DnaCharacterOrN> {
match kind.descendant {
TemplateSwitchDescendant::Query => {
get_char(&self.query, self.ranges.query_offset() + d - 1)
}
TemplateSwitchDescendant::Reference => {
get_char(&self.reference, self.ranges.reference_offset() + d - 1)
}
}
}
/// Offset-cost parameters for `kind`, precomputed once per secondary fill.
///
/// `compute_cost` charges `offset_costs(first_offset)` once at the template-switch entrance,
/// where `first_offset = a - r - o` for a jump from a pre_ts entrance at ancestor index `r`
/// into a secondary cell at logical ancestor index `a` (see [`Self::fill_jump_row`]). `o`
/// folds the ancestor range offset and the forward `-1` convention (see the `first_offset`
/// derivation in [`Self::retrace_ts`]). `pieces` decomposes the cost function into its finite
/// constant steps `(lo, hi, cost)`, which partition the offset support (one piece for the
/// common single-window config; several for a multi-step V-shape).
fn offset_model(&self, kind: TsKind) -> OffsetModel {
let costs = self.costs.offset_costs(kind.descendant, kind.ancestor);
let anc_offset = match kind.ancestor {
TemplateSwitchAncestor::Reference => self.ranges.reference_offset(),
TemplateSwitchAncestor::Query => self.ranges.query_offset(),
} as isize;
let o = anc_offset
+ match kind.direction {
TemplateSwitchDirection::Forward => 1,
TemplateSwitchDirection::Reverse => 0,
};
// A breakpoint `(input, cost)` sets the cost for offsets `[input, next_input - 1]` (the
// last extends to +infinity). Keep the finite ones as pieces; the unbounded sides are
// marked `isize::MIN` / `isize::MAX`.
let max = U64Cost::from_primitive(u64::MAX);
let breakpoints: Vec<(isize, U64Cost)> = costs.clone().into();
let pieces = breakpoints
.iter()
.enumerate()
.filter(|&(_, &(_, cost))| cost < max)
.map(|(i, &(lo, cost))| {
let hi = breakpoints
.get(i + 1)
.map_or(isize::MAX, |&(next, _)| next - 1);
(lo, hi, cost)
})
.collect();
OffsetModel { o, pieces }
}
/// Fill `jump_score` / `jump_r` for descendant column `desc_idx` (logical descendant
/// `desc_idx + 1`): for each secondary ancestor index `a` in `1..=a_len`, the cheapest pre_ts
/// entrance score *including* the offset cost, and the ancestor-axis index `r` achieving it.
/// `jump_score[a] == MAX` means no entrance lands in support for that `a`.
///
/// `e[r]` holds the pre_ts entrance scores along the ancestor axis at `desc_idx`. For each
/// offset-cost piece, the feasible `r` form the fixed-width window `[a - o - hi, a - o - lo]`
/// which slides right by one as `a` grows, so a monotone-index deque yields the windowed
/// minimum of `e` in amortised `O(1)` per `a`. Tie-break is lowest `r` (deque pops the back
/// on strict `>` only; the cross-piece merge prefers the smaller `r`), matching a plain
/// ascending scan over the support.
fn fill_jump_row(
model: &OffsetModel,
e: &[U64Cost],
a_len: usize,
jump_score: &mut [U64Cost],
jump_r: &mut [usize],
deque: &mut VecDeque<usize>,
) {
let max = U64Cost::from_primitive(u64::MAX);
for s in jump_score[..=a_len].iter_mut() {
*s = max;
}
if e.is_empty() {
return;
}
let last = (e.len() - 1) as isize;
for &(lo_p, hi_p, c_p) in &model.pieces {
deque.clear();
let mut r_added: isize = 0; // next ancestor index not yet pushed
for a in 1..=a_len {
let a_i = a as isize;
// first_offset t = a - r - o must lie in [lo_p, hi_p]:
// t <= hi_p => r >= a - o - hi_p (lower r bound)
// t >= lo_p => r <= a - o - lo_p (upper r bound)
let hi_r = if lo_p == isize::MIN {
last
} else {
a_i.saturating_sub(model.o).saturating_sub(lo_p).min(last)
};
let lo_r = if hi_p == isize::MAX {
0
} else {
a_i.saturating_sub(model.o).saturating_sub(hi_p).max(0)
};
if hi_r < 0 || lo_r > hi_r {
continue; // window empty for this `a` (it widens as `a` grows)
}
// Push indices entering the window from the right; pop dominated tail entries.
// Strict `>` keeps the lower index on value ties.
while r_added <= hi_r {
let ri = r_added as usize;
while deque.back().is_some_and(|&b| e[b] > e[ri]) {
deque.pop_back();
}
deque.push_back(ri);
r_added += 1;
}
// Drop indices that fell off the left.
let lo_r_u = lo_r as usize;
while deque.front().is_some_and(|&f| f < lo_r_u) {
deque.pop_front();
}
if let Some(&r) = deque.front() {
let total = e[r] + c_p;
if total < jump_score[a] || (total == jump_score[a] && r < jump_r[a]) {
jump_score[a] = total;
jump_r[a] = r;
}
}
}
}
}
/// Generic secondary fill, shared by all four template-switch matrices. Works in logical
/// `(ancestor a, descendant d)` coordinates and maps to physical cells via `kind.phys`.
fn fill_ts(&self, kind: TsKind, mat: &mut Mat<TsCell>) -> anyhow::Result<()> {
let (a_len, d_len) = self.ts_extents(kind);
let step = kind.ancestor_step();
let base = self
.costs
.base_cost
.get(kind.descendant, kind.ancestor, kind.direction);
let cost_table = match kind.direction {
TemplateSwitchDirection::Forward => &self.costs.secondary_forward_edit_costs,
TemplateSwitchDirection::Reverse => &self.costs.secondary_reverse_edit_costs,
};
let offset_model = self.offset_model(kind);
// Ancestor-axis length of pre_ts (rows for QR, cols for RQ), and buffers reused across
// descendant columns: the entrance-score row `e`, the per-`a` best entrance, and the
// sliding-window deque.
let max = U64Cost::from_primitive(u64::MAX);
let l = if kind.query_descendant() {
self.pre_ts.m()
} else {
self.pre_ts.n()
};
let mut e = vec![max; l];
let mut jump_score = vec![max; a_len + 1];
let mut jump_r = vec![0usize; a_len + 1];
let mut deque: VecDeque<usize> = VecDeque::new();
for d in 1..=d_len {
// Load the pre_ts entrance scores along the ancestor axis at descendant column
// `d - 1`, then precompute the cheapest offset-charged entrance for every `a`.
for (r, slot) in e.iter_mut().enumerate() {
let (row, col) = if kind.query_descendant() {
(r, d - 1)
} else {
(d - 1, r)
};
*slot = self.pre_ts[(row, col)].score;
}
Self::fill_jump_row(
&offset_model,
&e,
a_len,
&mut jump_score,
&mut jump_r,
&mut deque,
);
for a in 1..=a_len {
let anc = self.ts_ancestor_char(kind, a)?;
let dsc = self.ts_descendant_char(kind, d)?;
let matches = anc == dsc;
let mm_cost = cost_table.match_or_substitution_cost(anc, dsc);
let (row, col) = kind.phys(a, d);
// Mis-/Match: extend the secondary from its predecessor `(a - step, d - 1)`.
// Only a *set* predecessor is a real secondary alignment to extend; an unset
// cell (no entrance reached it — e.g. its offset is out of support) must not be
// mistaken for a length-0 chain.
let pred_a = a as isize - step;
let pred = (d > 1 && pred_a >= 1 && pred_a as usize <= a_len)
.then(|| kind.phys(pred_a as usize, d - 1))
.filter(|&(prow, pcol)| mat[(prow, pcol)].cell.ty.is_some());
if let Some((prow, pcol)) = pred {
let pred_score = mat[(prow, pcol)].cell.score;
let pred_len = mat[(prow, pcol)].ts_len;
let score = pred_score + mm_cost;
let cur = &mut mat[(row, col)];
if cur.cell.ty.is_none() || score < cur.cell.score {
cur.cell.ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
cur.cell.prev = None;
cur.cell.score = score;
cur.ts_len = pred_len + 1;
}
}
// Jump from pre_ts (template-switch entrance). The cheapest offset-charged
// entrance for this `a` was precomputed in `jump_score` / `jump_r`; the offset
// cost makes it depend on `a`, hence the per-column precompute.
if jump_score[a] < max {
let r = jump_r[a];
let entrance = if kind.query_descendant() {
(r, d - 1)
} else {
(d - 1, r)
};
let score = jump_score[a] + mm_cost + base;
let cur = &mut mat[(row, col)];
if cur.cell.ty.is_none() || score < cur.cell.score {
cur.cell.ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
cur.cell.prev = jump_to(WhichMatrix::PreTs, entrance.0, entrance.1);
cur.cell.score = score;
cur.ts_len = 1;
}
}
}
}
Ok(())
}
/// Cheapest secondary cell for `kind` at descendant index `desc_idx` (the descendant axis
/// fixed, scanning the ancestor axis). Returns the physical `(row, col)`, score, and inner
/// TS length (ties broken by longer TS). A skipped (empty) matrix reports a MAX score.
fn ts_best_at_descendant(&self, kind: TsKind, desc_idx: usize) -> (usize, usize, U64Cost, u16) {
let mat = self.ts_mat(kind.which);
let mut best_score = U64Cost::from_primitive(u64::MAX);
let mut best_len = 0u16;
let mut best = (0usize, 0usize);
let mut consider = |row: usize, col: usize, score: U64Cost, len: u16| {
if score < best_score || (score == best_score && len > best_len) {
best_score = score;
best_len = len;
best = (row, col);
}
};
if kind.query_descendant() {
// descendant = column: scan the ancestor rows.
for row in 1..mat.m() {
let c = &mat[(row, desc_idx)];
if c.cell.ty.is_some() {
consider(row, desc_idx, c.cell.score, c.ts_len);
}
}
} else {
// descendant = row: scan the ancestor columns.
for col in 1..mat.n() {
let c = &mat[(desc_idx, col)];
if c.cell.ty.is_some() {
consider(desc_idx, col, c.cell.score, c.ts_len);
}
}
}
(best.0, best.1, best_score, best_len)
}
/// Cheapest *complete* secondary cell for `kind`: a complete TS reaches the last descendant
/// index. A skipped (empty) matrix reports a MAX score.
fn ts_completion(&self, kind: TsKind) -> (usize, usize, U64Cost, u16) {
let mat = self.ts_mat(kind.which);
let last_descendant = if kind.query_descendant() {
mat.n().saturating_sub(1)
} else {
mat.m().saturating_sub(1)
};
self.ts_best_at_descendant(kind, last_descendant)
}
fn fill_post_ts(&mut self) -> anyhow::Result<()> {
let max_cost = U64Cost::from_primitive(u64::MAX);
// For each secondary matrix, precompute the cheapest jump source into post_ts, indexed
// by the post_ts descendant axis (column for QR, row for RQ). post_ts resumes one step
// past the template switch, so post_ts descendant index `k` reads the secondary's
// completion at descendant index `k - 1`; `k == 0` is a boundary and reports MAX.
// Skipped (empty) matrices report MAX everywhere and are never selected.
let sources: [Vec<(usize, usize, U64Cost)>; ALL_TS_KINDS.len()] =
std::array::from_fn(|idx| {
let kind = ALL_TS_KINDS[idx];
let len = if kind.query_descendant() {
self.post_ts.n()
} else {
self.post_ts.m()
};
(0..len)
.map(|k| {
if k >= 1 {
let (r, c, s, _) = self.ts_best_at_descendant(kind, k - 1);
(r, c, s)
} else {
(0, 0, max_cost)
}
})
.collect()
});
for j in 2..self.post_ts.n() {
for i in 1..self.post_ts.m() {
let rc = get_char(&self.reference, self.ranges.reference_offset() + i - 1)?;
let qc = get_char(&self.query, self.ranges.query_offset() + j - 1)?;
let matches = rc == qc;
let mm_cost = self
.costs
.primary_edit_costs
.match_or_substitution_cost(rc, qc);
// Template-switch exits: jump from each active secondary matrix. The base cost
// was already paid at the entrance. Iterating in `ALL_TS_KINDS` order keeps the
// earlier matrix on equal-score ties (the update test is strict `<`).
for (idx, kind) in ALL_TS_KINDS.iter().enumerate() {
let k = if kind.query_descendant() { j } else { i };
let (src_row, src_col, src_score) = sources[idx][k];
if src_score != max_cost {
let score = src_score + mm_cost;
let cur = &mut self.post_ts[(i, j)];
if cur.ty.is_none() || score < cur.score {
cur.ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
cur.prev = jump_to(kind.which, src_row, src_col);
cur.score = score;
}
}
}
if j > 2 && i > 1 {
// Mis-/Match
let cost = mm_cost;
let score = self.post_ts[(i - 1, j - 1)].score + cost;
if self.post_ts[(i, j)].ty.is_none() || score < self.post_ts[(i, j)].score {
self.post_ts[(i, j)].score = self.post_ts[(i - 1, j - 1)].score + cost;
self.post_ts[(i, j)].ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
self.post_ts[(i, j)].prev = None;
}
}
if i > 1 {
// Deletion
let cost = self.costs.primary_edit_costs.gap_extend_cost(rc);
let score = self.post_ts[(i - 1, j)].score + cost;
if self.post_ts[(i, j)].ty.is_none() || score < self.post_ts[(i, j)].score {
self.post_ts[(i, j)].score = score;
self.post_ts[(i, j)].ty = Some(Align::Deletion);
self.post_ts[(i, j)].prev = None;
}
}
if j > 2 {
// Insertion
let cost = self.costs.primary_edit_costs.gap_extend_cost(qc);
let score = self.post_ts[(i, j - 1)].score + cost;
if self.post_ts[(i, j)].ty.is_none() || score < self.post_ts[(i, j)].score {
self.post_ts[(i, j)].score = score;
self.post_ts[(i, j)].ty = Some(Align::Insertion);
self.post_ts[(i, j)].prev = None;
}
}
}
}
Ok(())
}
fn retrace_path(&self) -> anyhow::Result<AlignmentWithCost> {
use WhichMatrix::{PostTs, PreTs, TsQrf, TsQrr, TsRqf, TsRqr};
let best_result_pre_ts = self
.pre_ts
.fin()
.map_or(U64Cost::from_primitive(u64::MAX), |c| c.score);
let best_result_post_ts = self
.post_ts
.fin()
.filter(|c| c.ty.is_some())
.map_or(U64Cost::from_primitive(u64::MAX), |c| c.score);
// Each secondary matrix completes at the last descendant index; `ts_completion` scans
// the ancestor axis there for the cheapest cell. Collapse the four into one "complete
// TS" slot, lowest score winning and ties broken by longer inner TS, so the stable
// sort below can't let a shorter equal-score TS jump ahead of a longer one. Empty
// (skipped) matrices report a MAX score and are never selected.
let (best_complete_ts_score, best_complete_ts_mat, best_complete_ts_i, best_complete_ts_j) = {
let mut best_score = U64Cost::from_primitive(u64::MAX);
let mut best_len = 0u16;
let mut best_mat = PreTs;
let mut best_ij = (0usize, 0usize);
for kind in ALL_TS_KINDS {
let (ci, cj, score, len) = self.ts_completion(kind);
if score < best_score || (score == best_score && len > best_len) {
best_score = score;
best_len = len;
best_mat = kind.which;
best_ij = (ci, cj);
}
}
(best_score, best_mat, best_ij.0, best_ij.1)
};
let (score, mut curr_mat, mut curr_i, mut curr_j) = {
let mut indices = [
// Prefer no-TS, then complete TS (longer inner), then partial TS (shorter inner)
(
best_result_pre_ts,
PreTs,
self.pre_ts.m() - 1,
self.pre_ts.n() - 1,
),
(
best_complete_ts_score,
best_complete_ts_mat,
best_complete_ts_i,
best_complete_ts_j,
),
// saturating_sub: post_ts is empty when no TS direction is active, in which
// case its score is MAX and it is never selected.
(
best_result_post_ts,
PostTs,
self.post_ts.m().saturating_sub(1),
self.post_ts.n().saturating_sub(1),
),
];
indices.sort_by_key(|(res, ..)| *res);
indices[0]
};
let mut aln = Alignment::new();
while !(curr_mat == PreTs && curr_i == 0 && curr_j == 0) {
match curr_mat {
TsQrr | TsRqr | TsQrf | TsRqf => {
let kind = ts_kind_of(curr_mat);
self.retrace_ts(kind, &mut curr_mat, &mut curr_i, &mut curr_j, &mut aln)?;
}
_ => {
self.retrace_pre_or_post_ts(&mut curr_mat, &mut curr_i, &mut curr_j, &mut aln)?;
}
}
}
let actual_aln = aln.reverse();
Ok(AlignmentWithCost::new(actual_aln, score))
}
/// Given a completed secondary cell at `(new_i, new_j)` in matrix `mat`, walk back to its
/// entrance cell (the `ts_len == 1` jump) and return that entrance's stored pre_ts
/// back-pointer as `(row, col)`. The entrance offset depends on the ancestor-traversal
/// direction: reverse QRR moves up in rows / down in cols, reverse RQR the opposite,
/// forward QRF/RQF move down in both.
fn ts_entrance(
&self,
mat: WhichMatrix,
new_i: usize,
new_j: usize,
) -> anyhow::Result<(usize, usize)> {
use WhichMatrix::{TsQrf, TsQrr, TsRqf, TsRqr};
let ts = match mat {
TsQrr => &self.ts_qrr,
TsRqr => &self.ts_rqr,
TsQrf => &self.ts_qrf,
TsRqf => &self.ts_rqf,
_ => anyhow::bail!("ts_entrance called on non-secondary matrix"),
};
let len = ts[(new_i, new_j)].ts_len as usize;
let (entrance_i, entrance_j) = match mat {
TsQrr => (new_i + len - 1, new_j - len + 1),
TsRqr => (new_i - len + 1, new_j + len - 1),
TsQrf | TsRqf => (new_i - len + 1, new_j - len + 1),
_ => unreachable!(),
};
match ts[(entrance_i, entrance_j)].cell.prev {
Some((_, row, col)) => Ok((row as usize, col as usize)),
None => anyhow::bail!(
"We should have found the entrance jump in {mat:?} at i={entrance_i}, j={entrance_j}"
),
}
}
fn retrace_pre_or_post_ts(
&self,
curr_mat: &mut WhichMatrix,
curr_i: &mut usize,
curr_j: &mut usize,
aln: &mut Alignment<AlignmentType>,
) -> anyhow::Result<()> {
use WhichMatrix::{PostTs, PreTs, TsQrf, TsQrr, TsRqf, TsRqr};
let curr_mat_ptr = if *curr_mat == PreTs {
&self.pre_ts
} else {
&self.post_ts
};
let cell = &curr_mat_ptr[(*curr_i, *curr_j)];
match (cell.ty, cell.prev) {
// Local moves (no jump).
(Some(Align::Match), None) => {
aln.push(AlignmentType::PrimaryMatch);
*curr_i -= 1;
*curr_j -= 1;
}
(Some(Align::Substitution), None) => {
aln.push(AlignmentType::PrimarySubstitution);
*curr_i -= 1;
*curr_j -= 1;
}
(Some(Align::Insertion), None) => {
aln.push(AlignmentType::PrimaryInsertion);
*curr_j -= 1;
}
(Some(Align::Deletion), None) => {
aln.push(AlignmentType::PrimaryDeletion);
*curr_i -= 1;
}
// Template-switch exit: a (mis)match that jumped back from a secondary matrix.
(Some(ty @ (Align::Match | Align::Substitution)), Some((mat, new_i, new_j))) => {
let new_i = new_i as usize;
let new_j = new_j as usize;
aln.push(if ty == Align::Match {
AlignmentType::PrimaryMatch
} else {
AlignmentType::PrimarySubstitution
});
// Locate the secondary entrance cell (the `ts_len == 1` jump) and read its
// back-pointer into pre_ts. For QR* the anti-descendant is the reference (use
// the entrance row); for RQ* it is the query (use the entrance column).
let (entrance_row, entrance_col) = self.ts_entrance(mat, new_i, new_j)?;
let anti_descendant_gap = match mat {
TsQrr | TsQrf => isize::try_from(*curr_i)? - isize::try_from(entrance_row)? - 1,
TsRqr | TsRqf => isize::try_from(*curr_j)? - isize::try_from(entrance_col)? - 1,
PreTs | PostTs => anyhow::bail!("Implementation error"),
};
aln.push(AlignmentType::TemplateSwitchExit {
anti_descendant_gap,
});
*curr_mat = mat;
*curr_i = new_i;
*curr_j = new_j;
}
_ => anyhow::bail!("Implementation error"),
}
Ok(())
}
/// Generic secondary retrace, shared by all four template-switch matrices. Walks one step
/// of the secondary alignment, mapping logical `(ancestor, descendant)` moves back to
/// physical `(curr_i, curr_j)` via `kind`.
fn retrace_ts(
&self,
kind: TsKind,
curr_mat: &mut WhichMatrix,
curr_i: &mut usize,
curr_j: &mut usize,
aln: &mut Alignment<AlignmentType>,
) -> anyhow::Result<()> {
use WhichMatrix::PreTs;
let cell = &self.ts_mat(kind.which)[(*curr_i, *curr_j)];
if aln.inner_mut().is_empty() {
// Starting the retrace inside this matrix: emit the exit (end of the alignment).
// anti_descendant_gap = SP4 - SP1 on the ancestor (== anti-descendant), where SP4
// is the end of the ancestor's cluster range.
let (entrance_row, entrance_col) = self.ts_entrance(kind.which, *curr_i, *curr_j)?;
let (anc_range_len, anc_entrance) = if kind.query_descendant() {
(self.ranges.reference_range().len(), entrance_row)
} else {
(self.ranges.query_range().len(), entrance_col)
};
let gap = isize::try_from(anc_range_len + 1)? - isize::try_from(anc_entrance)? - 1;
aln.push(AlignmentType::TemplateSwitchExit {
anti_descendant_gap: gap,
});
}
let step = kind.ancestor_step();
match (cell.cell.ty, cell.cell.prev) {
// Secondary (mis)match extending within this matrix: step to its predecessor,
// logical `(a - step, d - 1)`.
(Some(ty @ (Align::Match | Align::Substitution)), None) => {
aln.push(if ty == Align::Match {
AlignmentType::SecondaryMatch
} else {
AlignmentType::SecondarySubstitution
});
let (a, d) = kind.unphys(*curr_i, *curr_j);
let pred_a = (a as isize - step) as usize;
let (nr, nc) = kind.phys(pred_a, d - 1);
*curr_i = nr;
*curr_j = nc;
}
// Template-switch entrance: jump back into pre_ts.
(Some(ty @ (Align::Match | Align::Substitution)), Some((PreTs, new_i, new_j))) => {
let new_i = new_i as usize;
let new_j = new_j as usize;
aln.push(if ty == Align::Match {
AlignmentType::SecondaryMatch
} else {
AlignmentType::SecondarySubstitution
});
// The ancestor absolute index is this matrix's ancestor axis at the current cell.
let (prim_entrance_asc, prim_offset, curr_anc_abs) = if kind.query_descendant() {
(new_i, self.ranges.reference_offset(), *curr_i)
} else {
(new_j, self.ranges.query_offset(), *curr_j)
};
// `first_offset` positions the secondary ancestor axis relative to the primary
// ancestor index at the entrance: `compute_cost` sets
// `ancestor_index = primary_ancestor + first_offset`. The first secondary char it
// then reads is `ancestor[ancestor_index]` for a forward TS but
// `ancestor[ancestor_index - 1].complement()` for a reverse TS, whereas our
// `ts_ancestor_char(a)` always maps logical `a` to `ancestor[a - 1]`. So the
// forward case needs an extra `-1` to line the two conventions up.
let mut first_offset = isize::try_from(curr_anc_abs)?
- isize::try_from(prim_entrance_asc)?
- isize::try_from(prim_offset)?;
if kind.direction == TemplateSwitchDirection::Forward {
first_offset -= 1;
}
aln.push(AlignmentType::TemplateSwitchEntrance {
first_offset,
equal_cost_range: EqualCostRange {
min_start: 0,
max_start: 0,
min_end: 0,
max_end: 0,
},
descendant: kind.descendant,
ancestor: kind.ancestor,
direction: kind.direction,
});
*curr_mat = PreTs;
*curr_i = new_i;
*curr_j = new_j;
}
_ => anyhow::bail!("Implementation error"),
}
Ok(())
}
}
/// Precomputed offset-cost parameters for one secondary fill (see [`AlignerInstance::offset_model`]).
struct OffsetModel {
/// Folds the ancestor range offset and the forward `-1` convention into `first_offset`.
o: isize,
/// Finite constant steps `(lo, hi, cost)` of the offset cost, partitioning its support.
/// `lo == isize::MIN` / `hi == isize::MAX` denote an unbounded side. Empty ⇒ no entrance
/// is ever affordable (the offset cost is everywhere infinite).
pieces: Vec<(isize, isize, U64Cost)>,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum WhichMatrix {
PreTs,
TsQrr,
TsRqr,
TsQrf,
TsRqf,
PostTs,
}
/// Describes one of the four template-switch matrices along the three axes that distinguish
/// them: which sequence is the descendant, which is the ancestor (== anti-descendant here),
/// and the copy direction. A single generic fill/retrace is driven by this descriptor.
///
/// Physical matrix layout is **not** normalised: a query-descendant (QR) matrix is laid out
/// `[ancestor(reference) rows] x [descendant(query) cols]`, while a reference-descendant (RQ)
/// matrix is the transpose `[descendant(reference) rows] x [ancestor(query) cols]`. The
/// `phys` mapping hides this so the generic code can work in logical `(ancestor, descendant)`
/// coordinates.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
struct TsKind {
which: WhichMatrix,
descendant: TemplateSwitchDescendant,
ancestor: TemplateSwitchAncestor,
direction: TemplateSwitchDirection,
}
const ALL_TS_KINDS: [TsKind; 4] = [
TsKind {
which: WhichMatrix::TsQrr,
descendant: TemplateSwitchDescendant::Query,
ancestor: TemplateSwitchAncestor::Reference,
direction: TemplateSwitchDirection::Reverse,
},
TsKind {
which: WhichMatrix::TsRqr,
descendant: TemplateSwitchDescendant::Reference,
ancestor: TemplateSwitchAncestor::Query,
direction: TemplateSwitchDirection::Reverse,
},
TsKind {
which: WhichMatrix::TsQrf,
descendant: TemplateSwitchDescendant::Query,
ancestor: TemplateSwitchAncestor::Reference,
direction: TemplateSwitchDirection::Forward,
},
TsKind {
which: WhichMatrix::TsRqf,
descendant: TemplateSwitchDescendant::Reference,
ancestor: TemplateSwitchAncestor::Query,
direction: TemplateSwitchDirection::Forward,
},
];
/// Recover the [`TsKind`] descriptor for one of the four secondary matrices.
fn ts_kind_of(which: WhichMatrix) -> TsKind {
ALL_TS_KINDS
.into_iter()
.find(|k| k.which == which)
.expect("ts_kind_of called on a non-secondary matrix")
}
impl TsKind {
/// True when the descendant is the query (the QR matrices), which are laid out
/// `ancestor(row) x descendant(col)`. The RQ matrices are the transpose.
const fn query_descendant(self) -> bool {
matches!(self.descendant, TemplateSwitchDescendant::Query)
}
/// Map logical `(ancestor, descendant)` indices to physical `(row, col)`.
const fn phys(self, ancestor: usize, descendant: usize) -> (usize, usize) {
if self.query_descendant() {
(ancestor, descendant)
} else {
(descendant, ancestor)
}
}
/// Inverse of [`Self::phys`]: physical `(row, col)` back to logical `(ancestor, descendant)`.
const fn unphys(self, row: usize, col: usize) -> (usize, usize) {
if self.query_descendant() {
(row, col)
} else {
(col, row)
}
}
/// Ancestor index step per descendant step: `+1` forward, `-1` reverse. The mis/match
/// predecessor of logical cell `(a, d)` is `(a - step, d - 1)`.
const fn ancestor_step(self) -> isize {
match self.direction {
TemplateSwitchDirection::Forward => 1,
TemplateSwitchDirection::Reverse => -1,
}
}
}
fn get_char(seq: &[u8], ix: usize) -> anyhow::Result<DnaCharacterOrN> {
let c = seq[ix];
DnaCharacterOrN::try_from(c)
.map_err(|()| anyhow!("Unknown character in query sequence: {}", &[c].as_bstr()))
}
/// Build a compact cross-matrix back-pointer (see [`Cell::prev`]). Matrix indices always
/// fit in `u32` (see the field docs), so the narrowing cast cannot lose information.
#[inline]
#[allow(clippy::cast_possible_truncation)]
fn jump_to(mat: WhichMatrix, i: usize, j: usize) -> Option<(WhichMatrix, u32, u32)> {
Some((mat, i as u32, j as u32))
}
#[derive(Debug, PartialEq, Eq, Clone, Copy)]
enum Align {
Match,
Substitution,
Deletion,
Insertion,
}
#[derive(Debug, PartialEq, Eq)]
struct Cell {
score: U64Cost,
ty: Option<Align>,
/// Predecessor of this cell.
/// `None` means the implicit local predecessor of the cell's own matrix (the diagonal
/// / row / column neighbour selected by `ty`). `Some((mat, i, j))` means this cell was
/// reached by a cross-matrix jump from `mat` at `(i, j)` (a template-switch entrance or
/// exit); `ty` still records whether that step is a (mis)match.
///
/// Indices are `u32` to keep this memory-critical struct compact: a matrix large enough
/// to overflow `u32` (billions of rows/cols) is far past what could ever be allocated.
/// The `Option` niches into `WhichMatrix`, so this field is 12 bytes, not 24.
prev: Option<(WhichMatrix, u32, u32)>,
}
impl Display for Cell {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let op = match self.ty {
Some(Align::Match) if self.prev.is_some() => "MaJ",
Some(Align::Substitution) if self.prev.is_some() => "SuJ",
Some(Align::Match) => "Mat",
Some(Align::Substitution) => "Sub",
Some(Align::Deletion) => "Del",
Some(Align::Insertion) => "Ins",
None => "???",
};
write!(f, "{:>3} ({op})", self.score.as_u64())
}
}
impl Default for Cell {
fn default() -> Self {
Self {
score: U64Cost::from_usize(0),
ty: None,
prev: None,
}
}
}
#[derive(Debug, Default, PartialEq, Eq)]
struct TsCell {
cell: Cell,
ts_len: u16,
}
impl Display for TsCell {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{} (len={:>3})", self.cell, self.ts_len)
}
}
#[derive(Clone)]
struct Mat<T> {
data: Vec<T>,
n: usize,
}
impl<T: Default> Mat<T> {
pub fn try_new(m: usize, n: usize) -> Result<Self, AlignmentFailure> {
assert!(n > 0, "width must be nonzero");
let mut data = Vec::new();
data.try_reserve_exact(m * n)
.map_err(|_| AlignmentFailure::oom())?;
data.resize_with(m * n, T::default);
Ok(Self { data, n })
}
}
impl<T> Mat<T> {
pub fn empty() -> Self {
Self {
data: Vec::new(),
n: 1,
}
}
pub const fn n(&self) -> usize {
self.n
}
pub const fn m(&self) -> usize {
self.data.len() / self.n
}
/// i = row index (< m)
/// j = col index (< n)
#[inline]
fn index(&self, i: usize, j: usize) -> usize {
debug_assert!(
i < self.m(),
"row {i} out of bounds, matrix has {} rows and {} columns",
self.m(),
self.n
);
debug_assert!(
j < self.n,
"column {j} out of bounds, matrix has {} rows and {} columns",
self.m(),
self.n
);
i * self.n + j
}
fn fin(&self) -> Option<&T> {
self.data.last()
}
}
impl<T> Index<(usize, usize)> for Mat<T> {
type Output = T;
fn index(&self, (i, j): (usize, usize)) -> &Self::Output {
let idx = self.index(i, j);
&self.data[idx]
}
}
impl<T> IndexMut<(usize, usize)> for Mat<T> {
fn index_mut(&mut self, (i, j): (usize, usize)) -> &mut Self::Output {
let idx = self.index(i, j);
&mut self.data[idx]
}
}
impl<T: fmt::Display> fmt::Debug for Mat<T> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let height = self.m();
writeln!(f, "[")?;
for i in 0..height {
write!(f, " [")?;
for j in 0..self.n {
write!(f, "{}", self[(i, j)])?;
if j + 1 < self.n {
write!(f, ", ")?;
}
}
writeln!(f, "],")?;
}
writeln!(f, "]")?;
Ok(())
}
}
#[cfg(test)]
mod tests {
use std::sync::Arc;
use lib_tsalign::a_star_aligner::alignment_geometry::AlignmentCoordinates;
use super::*;
fn has_ts(geom: &str, alignment: &TwitcherAlignmentWithStatistics) -> bool {
assert_eq!(geom.len(), 3);
let lowercase = geom.to_lowercase();
let mut geom = lowercase.chars();
let p = match geom.next().unwrap() {
'r' => TemplateSwitchDescendant::Reference,
'q' => TemplateSwitchDescendant::Query,
_ => panic!(),
};
let s = match geom.next().unwrap() {
'r' => TemplateSwitchAncestor::Reference,
'q' => TemplateSwitchAncestor::Query,
_ => panic!(),
};
let d = match geom.next().unwrap() {
'r' => TemplateSwitchDirection::Reverse,
'f' => TemplateSwitchDirection::Forward,
_ => panic!(),
};
for (_, ty) in alignment.alignment.alignment.iter_compact() {
match ty {
AlignmentType::TemplateSwitchEntrance {
descendant,
ancestor,
direction,
..
} if *descendant == p && *ancestor == s && *direction == d => return true,
_ => {}
}
}
false
}
fn cig(r: &[u8], q: &[u8], range: AlignmentRange) -> String {
let fpa = FourPointAligner::default();
let res = fpa.align(Arc::from(r), Arc::from(q), range).unwrap();
res.alignment.alignment.cigar()
}
/// Replays an fpa alignment over the sequences using the exact coordinate conventions of
/// `Alignment::compute_cost` (lib_tsalign) and asserts that every emitted (mis)match label is
/// consistent with the bases it lands on. This validates the TS geometry — `first_offset`,
/// `anti_descendant_gap`, and the descendant/ancestor/direction roles — independent of the
/// cost model (fpa deliberately omits length/offset/gap penalties, so costs may differ; the
/// *character alignment* must not).
fn assert_geometry(r: &[u8], q: &[u8], range: AlignmentRange) {
let (r_off, q_off) = (range.reference_offset(), range.query_offset());
let fpa = FourPointAligner::default();
let res = fpa.align(Arc::from(r), Arc::from(q), range).unwrap();
let cigar = res.alignment.alignment.cigar();
let chr = |seq: &[u8], i: usize| DnaCharacterOrN::try_from(seq[i]).unwrap();
let to_isize = |x: usize| isize::try_from(x).unwrap();
let mut ri = r_off;
let mut qi = q_off;
let mut di = 0usize; // descendant_index
let mut ai = 0usize; // ancestor_index
let mut desc = TemplateSwitchDescendant::Reference;
let mut anc = TemplateSwitchAncestor::Reference;
let mut dir = TemplateSwitchDirection::Forward;
for ty in res.alignment.alignment.iter_flat_cloned() {
match ty {
AlignmentType::PrimaryInsertion => qi += 1,
AlignmentType::PrimaryDeletion => ri += 1,
AlignmentType::PrimaryMatch | AlignmentType::PrimarySubstitution => {
let is_match = chr(r, ri) == chr(q, qi);
assert_eq!(
is_match,
ty == AlignmentType::PrimaryMatch,
"primary label/base inconsistency at r={ri} q={qi} in {cigar}"
);
ri += 1;
qi += 1;
}
AlignmentType::TemplateSwitchEntrance {
first_offset,
descendant,
ancestor,
direction,
..
} => {
desc = descendant;
anc = ancestor;
dir = direction;
di = match desc {
TemplateSwitchDescendant::Reference => ri,
TemplateSwitchDescendant::Query => qi,
};
let prim_anc = match anc {
TemplateSwitchAncestor::Reference => ri,
TemplateSwitchAncestor::Query => qi,
};
ai = usize::try_from(to_isize(prim_anc) + first_offset).unwrap();
}
AlignmentType::SecondaryMatch | AlignmentType::SecondarySubstitution => {
let dc = match desc {
TemplateSwitchDescendant::Reference => chr(r, di),
TemplateSwitchDescendant::Query => chr(q, di),
};
let ac = match (anc, dir) {
(TemplateSwitchAncestor::Reference, TemplateSwitchDirection::Forward) => {
chr(r, ai)
}
(TemplateSwitchAncestor::Reference, TemplateSwitchDirection::Reverse) => {
chr(r, ai - 1).complement()
}
(TemplateSwitchAncestor::Query, TemplateSwitchDirection::Forward) => {
chr(q, ai)
}
(TemplateSwitchAncestor::Query, TemplateSwitchDirection::Reverse) => {
chr(q, ai - 1).complement()
}
};
assert_eq!(
dc == ac,
ty == AlignmentType::SecondaryMatch,
"secondary label/base inconsistency at descendant={di} ancestor={ai} in {cigar}"
);
di += 1;
match dir {
TemplateSwitchDirection::Forward => ai += 1,
TemplateSwitchDirection::Reverse => ai -= 1,
}
}
AlignmentType::SecondaryInsertion => di += 1,
AlignmentType::SecondaryDeletion => match dir {
TemplateSwitchDirection::Forward => ai += 1,
TemplateSwitchDirection::Reverse => ai -= 1,
},
AlignmentType::TemplateSwitchExit {
anti_descendant_gap,
} => match desc {
TemplateSwitchDescendant::Reference => {
ri = di;
qi = usize::try_from(to_isize(qi) + anti_descendant_gap).unwrap();
}
TemplateSwitchDescendant::Query => {
qi = di;
ri = usize::try_from(to_isize(ri) + anti_descendant_gap).unwrap();
}
},
other => panic!("unexpected alignment type {other:?} in {cigar}"),
}
}
}
#[test]
fn geometry_qrr_rqr() {
// Reverse TS in both descendant orientations (the vcf-style cases).
assert_geometry(
b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
b"AAAAAAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAAA",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(15, 15),
AlignmentCoordinates::new(26, 25),
),
);
assert_geometry(
b"AAAAAAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAAA",
b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(15, 15),
AlignmentCoordinates::new(25, 26),
),
);
}
#[test]
fn geometry_qrf_rqf() {
const FLANK_L: &str = "GGGGGGGGGGGGGGG";
const SEG_A: &str = "TACGTTCAGGAC";
const FLANK_M: &str = "AAAAAAAAAAAAAAA";
const SEG_B: &str = "CTCTCTCTCTCT";
const FLANK_R: &str = "TTTTTTTTTTTTTTT";
let with_a = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_A}{FLANK_R}");
let with_b = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_B}{FLANK_R}");
// QRF: forward duplication in the query.
assert_geometry(
with_b.as_bytes(),
with_a.as_bytes(),
AlignmentRange::new_complete(with_b.len(), with_a.len()),
);
// RQF: forward duplication in the reference.
assert_geometry(
with_a.as_bytes(),
with_b.as_bytes(),
AlignmentRange::new_complete(with_a.len(), with_b.len()),
);
}
#[test]
fn geometry_large() {
let r = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTTTTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTCAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCCACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGAGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAATACTAGATAGAGAATTCATCCGCTAACCTTTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGCAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCACAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCCGTTCTAACATTCTGTAATTGGGTAGTACTGGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTATGGTCAAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
let q = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTATTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTTAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCTACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGCGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAGTACTAGATAGAAAGTTAGCGGATGAATTCTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACCATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGAAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCGCAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCTGTTCTAACATTCTGTAATTGGGTAGTACTGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTACGGTCAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
assert_geometry(r, q, AlignmentRange::new_complete(r.len(), q.len()));
assert_geometry(q, r, AlignmentRange::new_complete(q.len(), r.len()));
}
/// The DP matrices hold one `Cell` / `TsCell` per (row, col) and dominate memory, so guard
/// their size against accidental growth. Critically, the `prev` back-pointer must niche
/// into `WhichMatrix` (12 bytes, not 16) and use `u32` indices; widening those to `usize`
/// would push `prev` to 24 and `TsCell` to 48.
#[test]
fn cell_layout_is_compact() {
use std::mem::size_of;
assert_eq!(
size_of::<Option<(WhichMatrix, u32, u32)>>(),
12,
"prev back-pointer lost its niche optimization"
);
assert_eq!(size_of::<Cell>(), 24, "Cell grew");
assert_eq!(size_of::<TsCell>(), 32, "TsCell grew");
}
#[test]
fn fpa_basic() {
assert_eq!(
cig(b"ACGTA", b"AACGTTA", AlignmentRange::new_complete(5, 7)),
"1I3=1I2="
);
}
#[test]
fn fpa_basic_inv() {
assert_eq!(
cig(b"AACGTTA", b"ACGTA", AlignmentRange::new_complete(7, 5)),
"1D3=1D2="
);
}
#[test]
fn fpa_with_ranges() {
assert_eq!(
cig(
b"TTTTTACGTATTTTT",
b"TTTTTACGAAAAAAATTTTT",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(5, 5),
AlignmentCoordinates::new(10, 15),
),
),
"1=[TSRQR:[0,0]:[0,0]:2:4=:9]"
);
}
#[test]
fn fpa_with_ranges_inv() {
assert_eq!(
cig(
b"TTTTTACGAAAAAAATTTTT",
b"TTTTTACGTATTTTT",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(5, 5),
AlignmentCoordinates::new(15, 10),
),
),
"1=[TSQRR:[0,0]:[0,0]:2:4=:9]"
);
}
#[test]
fn fpa_with_cursed_ranges() {
assert_eq!(
cig(
b"TTTTTACGTATTTTT",
b"TTTTTACGAAAAAAATTTTT",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(0, 18),
AlignmentCoordinates::new(1, 20),
),
),
"1I1="
);
}
#[test]
fn fpa_with_cursed_ranges_inv() {
assert_eq!(
cig(
b"TTTTTACGAAAAAAATTTTT",
b"TTTTTACGTATTTTT",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(18, 0),
AlignmentCoordinates::new(20, 1),
),
),
"1D1="
);
}
#[test]
fn fpa_vcf() {
assert_eq!(
cig(
b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
b"AAAAAAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAAA",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(15, 15),
AlignmentCoordinates::new(26, 25),
),
),
"[TSQRR:[0,0]:[0,0]:-5:10=:11]"
);
}
#[test]
fn fpa_vcf_inv() {
assert_eq!(
cig(
b"AAAAAAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAAA",
b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(15, 15),
AlignmentCoordinates::new(25, 26),
),
),
"[TSRQR:[0,0]:[0,0]:-5:10=:11]"
);
}
#[test]
fn fpa_large() {
let r = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTTTTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTCAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCCACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGAGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAATACTAGATAGAGAATTCATCCGCTAACCTTTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGCAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCACAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCCGTTCTAACATTCTGTAATTGGGTAGTACTGGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTATGGTCAAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
let q = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTATTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTTAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCTACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGCGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAGTACTAGATAGAAAGTTAGCGGATGAATTCTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACCATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGAAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCGCAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCTGTTCTAACATTCTGTAATTGGGTAGTACTGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTACGGTCAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
let fpa = FourPointAligner::default();
let res = fpa
.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_complete(r.len(), q.len()),
)
.unwrap();
assert!(has_ts("qrr", &res));
assert_eq!(
res.alignment.alignment.cigar(),
"125=1X114=1X124=1X42=1X92=1X11=1D2=[TSQRR:[0,0]:[0,0]:13:23=:23]36=1I94=1I88=1X43=1X61=1X29=1D92=1X4=1D39="
);
}
#[test]
fn fpa_large_inv() {
let q = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTTTTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTCAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCCACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGAGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAATACTAGATAGAGAATTCATCCGCTAACCTTTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGCAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCACAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCCGTTCTAACATTCTGTAATTGGGTAGTACTGGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTATGGTCAAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
let r = b"TTTGTTTTCTTTTTAGGAAGCCATCATTCTTTAGAGGGCAATGACCAAACAGTACCAGCAGAAATTGAAGTACCAGCAGAAGGCTAAGAAGGTTAGGAGAAAAAGACATTTTGTATTTTACTGTTATTTTCTTTTTTTTTTTTTAGACGAAGTCTTGCTCTGTCACCAGGCTAGAGTGCAGTGGCACGATCTCGGCTCACTGCAACCTCTGACTCCCTGGTTCAAGTGATTCTCCTGCCTTAGCCTCCTAAGTAGCTGGGATTATAGGCACAGACCACCACATCCAGCTATTTTTTGTATTTTTGGTAGAGACCAGGGTTTCACCATGTTGGCCAGGATGGTCTCAATCTTTTGACCTCCTGATCTACCCACCTCGGCCTCCCAACATGCTGGGATTACAGATGTGGGCGCTTGGCCACCTCCTCTTGGGAGAAATGCACTGATTCTGGTTGCCACGTGGATTTATTTTGGGAGTGATATTCATCTAACTTCATGGAAATAGTACTAGATAGAAAGTTAGCGGATGAATTCTCTATCTGATGAGAGTTTTGGGCAAATCGAATACCAAGTTACCAAGTTTTGTTTTTTTCTCTGATGCAAAAAAACAATTTGCCAGCCAGTGAAAAACTCTCACAGCTCTGGATGTGAGTTTAGGATACTGGATTTCTACCATTCAATTTCTTACTACTTTTCTTGCACAGGGATCATGGCACAAGCTGCAGTTTCCACCCTGCCCATTGAAGATGAGGAGTCCATGGAAGATGAGGAGTCCGTTGAAGATGAGTCTGTTGAAGATGAGTCCGCAGAGAACAGGATGGTGGTGACATTGCTCATATCAGCTCTTGAGTCCATGGTGAGACCTTCTGTTCTAACATTCTGTAATTGGGTAGTACTGGGTGGTAGATAAGGTTGATTTGTTTTTGTAGAATTTATAATTTTATGATTTATAGTTCTAATGAGTAGATCTTTTTCTTGAATAGTAGTTACGGTCAAACACTTCTGACCAAATGTGCCATGTTGTCCAGCCTGG";
let fpa = FourPointAligner::default();
let res = fpa
.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_complete(r.len(), q.len()),
)
.unwrap();
assert!(has_ts("rqr", &res));
assert_eq!(
res.alignment.alignment.cigar(),
"125=1X114=1X124=1X42=1X92=1X11=1I2=[TSRQR:[0,0]:[0,0]:13:23=:23]36=1D94=1D88=1X43=1X61=1X29=1I92=1X4=1I39="
);
}
// Forward (same-strand) template switch: the query contains a forward copy of an
// earlier reference segment (a direct duplication) in place of the unrelated segB.
// A monotonic primary alignment cannot reuse segA, so the cheapest explanation is a
// forward TS (descendant=Query, ancestor=Reference).
#[test]
fn fpa_forward_qrf() {
const FLANK_L: &str = "GGGGGGGGGGGGGGG";
const SEG_A: &str = "TACGTTCAGGAC";
const FLANK_M: &str = "AAAAAAAAAAAAAAA";
const SEG_B: &str = "CTCTCTCTCTCT";
const FLANK_R: &str = "TTTTTTTTTTTTTTT";
let r = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_B}{FLANK_R}");
let q = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_A}{FLANK_R}");
let fpa = FourPointAligner::default();
let res = fpa
.align(
Arc::from(r.as_bytes()),
Arc::from(q.as_bytes()),
AlignmentRange::new_complete(r.len(), q.len()),
)
.unwrap();
assert!(has_ts("qrf", &res));
assert_eq!(
res.alignment.alignment.cigar(),
"42=[TSQRF:[0,0]:[0,0]:-27:12=:12]15="
);
}
const QRF_FLANK_L: &str = "GGGGGGGGGGGGGGG";
const QRF_SEG_A: &str = "TACGTTCAGGAC";
const QRF_FLANK_M: &str = "AAAAAAAAAAAAAAA";
const QRF_SEG_B: &str = "CTCTCTCTCTCT";
const QRF_FLANK_R: &str = "TTTTTTTTTTTTTTT";
/// Build the `fpa_forward_qrf` scenario: a forward duplication of `SEG_A` in the query,
/// explainable only by a QRF template switch at `first_offset == -27` (the single reference
/// copy of `SEG_A`).
fn qrf_seqs() -> (String, String) {
let r = format!("{QRF_FLANK_L}{QRF_SEG_A}{QRF_FLANK_M}{QRF_SEG_B}{QRF_FLANK_R}");
let q = format!("{QRF_FLANK_L}{QRF_SEG_A}{QRF_FLANK_M}{QRF_SEG_A}{QRF_FLANK_R}");
(r, q)
}
/// Align with a custom `rq_qr` offset cost function (used by the QR/RQ template switches),
/// everything else the default config.
fn align_with_rq_qr_offset(
r: &[u8],
q: &[u8],
range: AlignmentRange,
rq_qr: Vec<(isize, U64Cost)>,
) -> TwitcherAlignmentWithStatistics {
use lib_tsalign::costs::cost_function::CostFunction;
let config = TemplateSwitchConfig::<DnaAlphabetOrN, U64Cost> {
rq_qr_offset_costs: CostFunction::try_from(rq_qr).unwrap(),
..Default::default()
};
let fpa = FourPointAligner::new(config, false);
fpa.align(Arc::from(r), Arc::from(q), range).unwrap()
}
// An offset cost of infinity outside a window must forbid template-switch entrances whose
// `first_offset` lands outside it. The only QRF explanation here sits at offset -27, so a
// window of [0, 100] (excluding -27) must suppress the QRF template switch entirely.
#[test]
fn offset_cost_forbids_out_of_window_qrf() {
let max = U64Cost::from_primitive(u64::MAX);
let zero = U64Cost::from_primitive(0);
let (r, q) = qrf_seqs();
let res = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![(isize::MIN, max), (0, zero), (101, max)],
);
assert!(!has_ts("qrf", &res), "offset -27 excluded yet QRF emitted");
}
// The complementary case: a window of exactly {-27} keeps the QRF template switch. Guards
// against over-restricting (the in-window entrance must still be reachable).
#[test]
fn offset_cost_keeps_in_window_qrf() {
let max = U64Cost::from_primitive(u64::MAX);
let zero = U64Cost::from_primitive(0);
let (r, q) = qrf_seqs();
let res = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![(isize::MIN, max), (-27, zero), (-26, max)],
);
assert!(has_ts("qrf", &res), "offset -27 allowed yet no QRF");
}
// A constant (flat) offset cost is V-shaped, does not change which entrance is optimal, and
// must be charged exactly once at the template-switch entrance. So raising the constant from
// 0 to 5 raises the total alignment cost by exactly 5, geometry unchanged.
#[test]
fn constant_offset_cost_adds_to_score_once() {
let zero = U64Cost::from_primitive(0);
let five = U64Cost::from_primitive(5);
let (r, q) = qrf_seqs();
let base = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![(isize::MIN, zero)],
);
let plus = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![(isize::MIN, five)],
);
assert!(has_ts("qrf", &base));
assert!(has_ts("qrf", &plus));
assert_eq!(
plus.alignment.cost,
base.alignment.cost + five,
"constant offset cost not charged exactly once"
);
}
// A multi-level V-shaped offset cost (several distinct finite cost steps) must charge the
// exact step that `first_offset` lands in. The only QRF source sits at offset -27, which the
// cost function below maps to cost 4 (its [-50,-11] step). Exercises the k>1 piece handling.
#[test]
fn multi_level_offset_cost_charges_landing_step() {
let max = U64Cost::from_primitive(u64::MAX);
let zero = U64Cost::from_primitive(0);
let four = U64Cost::from_primitive(4);
let (r, q) = qrf_seqs();
// V-shape: MAX | [-50,-11]=4 | [-10,10]=0 | [11,100]=4 | MAX. offset -27 -> 4.
let base = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![(isize::MIN, zero)],
);
let stepped = align_with_rq_qr_offset(
r.as_bytes(),
q.as_bytes(),
AlignmentRange::new_complete(r.len(), q.len()),
vec![
(isize::MIN, max),
(-50, four),
(-10, zero),
(11, four),
(101, max),
],
);
assert!(has_ts("qrf", &base));
assert!(has_ts("qrf", &stepped));
assert_eq!(
stepped.alignment.cost,
base.alignment.cost + four,
"offset -27 not charged its V-shape step (4)"
);
}
// Same as above but the forward duplication lives in the reference, copied from the
// query (descendant=Reference, ancestor=Query).
#[test]
fn fpa_forward_rqf() {
const FLANK_L: &str = "GGGGGGGGGGGGGGG";
const SEG_A: &str = "TACGTTCAGGAC";
const FLANK_M: &str = "AAAAAAAAAAAAAAA";
const SEG_B: &str = "CTCTCTCTCTCT";
const FLANK_R: &str = "TTTTTTTTTTTTTTT";
let q = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_B}{FLANK_R}");
let r = format!("{FLANK_L}{SEG_A}{FLANK_M}{SEG_A}{FLANK_R}");
let fpa = FourPointAligner::default();
let res = fpa
.align(
Arc::from(r.as_bytes()),
Arc::from(q.as_bytes()),
AlignmentRange::new_complete(r.len(), q.len()),
)
.unwrap();
assert!(has_ts("rqf", &res));
assert_eq!(
res.alignment.alignment.cigar(),
"42=[TSRQF:[0,0]:[0,0]:-27:12=:12]15="
);
}
}