use std::{
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, TemplateSwitchPrimary},
},
config::TemplateSwitchConfig,
costs::U64Cost,
};
use serde::{Deserialize, Serialize};
use tracing::instrument;
use crate::common::{
ImmutableSequence,
aligner::result::{
AlignmentFailure, TwitcherAlignment, TwitcherAlignmentCase, TwitcherAlignmentResult,
},
};
#[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,
mat1: Mat<Cell>,
mat2: Mat<RevCell>,
mat3: Mat<Cell>,
}
#[allow(dead_code)]
#[derive(Debug, Clone, Deserialize, Serialize)]
pub struct FpaAlignmentStatistics {
pub duration: Duration,
pub ranges: AlignmentRange,
}
impl FourPointAligner {
pub 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> {
if ranges.reference_range().is_empty() {
return Err(AlignmentFailure::error(
"Cannot align if the reference range is empty",
));
}
if ranges.query_range().len() < 2 {
return Err(AlignmentFailure::error(
"Cannot align if the length of the query range is less than 2",
));
}
let start_time = Instant::now();
// First matrix only aligns within range
let mat1: Mat<Cell> = Mat::try_new(
ranges.reference_range().len() + 1,
ranges.query_range().len() + 1,
)?;
// Second Matrix aligns in the whole reference sequence, but only the relevant path of the query sequence. This is the big one in our case.
// TODO since we only consider matches, surely this can be done in linear space at least... :D
let mat2: Mat<RevCell> = Mat::try_new(reference.len() + 1, ranges.query_range().len() + 1)?;
// Third matrix again only aligns within range
let mat3: Mat<Cell> = Mat::try_new(
ranges.reference_range().len() + 1,
ranges.query_range().len() + 1,
)?;
Ok(Self {
start_time,
reference,
query,
ranges,
costs,
no_ts,
mat1,
mat2,
mat3,
})
}
fn run(mut self) -> Result<TwitcherAlignment, AlignmentFailure> {
self.fill_mat1()?;
if !self.no_ts {
self.fill_mat2()?;
self.fill_mat3()?;
}
let ts_alignment = self.retrace_path()?;
let statistics = FpaAlignmentStatistics {
duration: self.start_time.elapsed(),
ranges: self.ranges,
};
Ok(TwitcherAlignment {
result: ts_alignment,
stats: statistics.into(),
})
}
fn fill_mat1(&mut self) -> anyhow::Result<()> {
for i in 0..self.mat1.m() {
for j in 0..self.mat1.n() {
let rc = (i > 0)
.then(|| get_char(&self.reference, self.ranges.reference_offset() + i - 1))
.transpose()?;
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.mat1[(i - 1, j - 1)].score + cost;
if self.mat1[(i, j)].ty.is_none() || score < self.mat1[(i, j)].score {
self.mat1[(i, j)].score = self.mat1[(i - 1, j - 1)].score + cost;
self.mat1[(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.mat1[(i - 1, j)].score + cost;
if self.mat1[(i, j)].ty.is_none() || score < self.mat1[(i, j)].score {
self.mat1[(i, j)].score = score;
self.mat1[(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.mat1[(i, j - 1)].score + cost;
if self.mat1[(i, j)].ty.is_none() || score < self.mat1[(i, j)].score {
self.mat1[(i, j)].score = score;
self.mat1[(i, j)].ty = Some(Align::Insertion);
}
}
}
}
Ok(())
}
fn fill_mat2(&mut self) -> anyhow::Result<()> {
for j in 1..self.mat2.n() {
let mut best_score = U64Cost::from_usize(usize::MAX);
let mut best_i = 0;
for i in 0..self.mat1.m() {
if self.mat1[(i, j - 1)].score < best_score {
best_score = self.mat1[(i, j - 1)].score;
best_i = i;
}
}
for i in 1..self.mat2.m() {
// Use indices directly as this matrix has the whole sequence
let rc = get_char(&self.reference, i - 1)?;
let rc_rc = rc.complement();
let qc = get_char(&self.query, self.ranges.query_offset() + j - 1)?;
let matches = rc_rc == qc;
let mm_cost = self
.costs
.secondary_reverse_edit_costs
.match_or_substitution_cost(rc_rc, qc);
if j > 1 && i < self.mat2.m() - 1 {
// Mis-/Match
let cost = mm_cost;
let score = self.mat2[(i + 1, j - 1)].cell.score + cost;
if self.mat2[(i, j)].cell.ty.is_none() || score < self.mat2[(i, j)].cell.score {
self.mat2[(i, j)].cell.ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
self.mat2[(i, j)].cell.score = score;
self.mat2[(i, j)].ts_len = self.mat2[(i + 1, j - 1)].ts_len + 1;
}
}
{
// Jump
let cost = mm_cost + self.costs.base_cost.qrr;
let score = self.mat1[(best_i, j - 1)].score + cost;
if self.mat2[(i, j)].cell.ty.is_none() || score < self.mat2[(i, j)].cell.score {
self.mat2[(i, j)].cell.ty = Some(if matches {
Align::MatchJump(best_i, j - 1)
} else {
Align::SubstitutionJump(best_i, j - 1)
});
self.mat2[(i, j)].cell.score = score;
self.mat2[(i, j)].ts_len = 1;
}
}
}
}
Ok(())
}
fn fill_mat3(&mut self) -> anyhow::Result<()> {
for j in 2..self.mat3.n() {
let mut best_score = U64Cost::from_usize(usize::MAX);
let mut best_i = 0;
// Note we start iterating at 1
for i in 1..self.mat2.m() {
if self.mat2[(i, j - 1)].cell.score < best_score {
best_score = self.mat2[(i, j - 1)].cell.score;
best_i = i;
}
}
for i in 1..self.mat3.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);
if j > 2 && i > 1 {
// Mis-/Match
let cost = mm_cost;
let score = self.mat3[(i - 1, j - 1)].score + cost;
if self.mat3[(i, j)].ty.is_none() || score < self.mat3[(i, j)].score {
self.mat3[(i, j)].score = self.mat3[(i - 1, j - 1)].score + cost;
self.mat3[(i, j)].ty = Some(if matches {
Align::Match
} else {
Align::Substitution
});
}
}
if i > 1 {
// Deletion
let cost = self.costs.primary_edit_costs.gap_extend_cost(rc);
let score = self.mat3[(i - 1, j)].score + cost;
if self.mat3[(i, j)].ty.is_none() || score < self.mat3[(i, j)].score {
self.mat3[(i, j)].score = score;
self.mat3[(i, j)].ty = Some(Align::Deletion);
}
}
if j > 2 {
// Insertion
let cost = self.costs.primary_edit_costs.gap_extend_cost(qc);
let score = self.mat3[(i, j - 1)].score + cost;
if self.mat3[(i, j)].ty.is_none() || score < self.mat3[(i, j)].score {
self.mat3[(i, j)].score = score;
self.mat3[(i, j)].ty = Some(Align::Insertion);
}
}
{
// Jump
let cost = mm_cost; // Because the jump cost has been added for the first jump already, we do not add it again
let score = self.mat2[(best_i, j - 1)].cell.score + cost;
if self.mat3[(i, j)].ty.is_none() || score < self.mat3[(i, j)].score {
self.mat3[(i, j)].ty = Some(if matches {
Align::MatchJump(best_i, j - 1)
} else {
Align::SubstitutionJump(best_i, j - 1)
});
self.mat3[(i, j)].score = score;
}
}
}
}
Ok(())
}
fn retrace_path(&mut self) -> anyhow::Result<TwitcherAlignmentCase> {
use CurrentMatrix::{One, Three, Two};
let best_result_mat1 = self.mat1[(self.mat1.m() - 1, self.mat1.n() - 1)].score;
let (best_result_mat2_i, best_result_mat2) = {
let mut best_i = 0;
let mut best_score: U64Cost = U64Cost::from_primitive(u64::MAX);
for i in 1..self.mat2.m() {
let cell = &self.mat2[(i, self.mat2.n() - 1)];
if cell.cell.score < best_score {
best_i = i;
best_score = cell.cell.score;
}
}
(best_i, best_score)
};
let best_result_mat3 = self.mat3[(self.mat3.m() - 1, self.mat3.n() - 1)].score;
let has_ts = best_result_mat3 < best_result_mat1 || best_result_mat2 < best_result_mat1;
let (_score, mut curr_mat, mut curr_i, mut curr_j) = {
let mut indices = [
// The order matters here: We prefer no ts, then a complete ts, and least prefer a partial ts (on equal score)
(best_result_mat1, One, self.mat1.m() - 1, self.mat1.n() - 1),
(best_result_mat2, Two, best_result_mat2_i, self.mat2.n() - 1),
(
best_result_mat3,
Three,
self.mat3.m() - 1,
self.mat3.n() - 1,
),
];
indices.sort_by_key(|(res, ..)| *res);
indices[0]
};
let mut aln = Alignment::new();
while !(curr_mat == One && curr_i == 0 && curr_j == 0) {
if curr_mat == Two {
self.retrace_mat2(&mut curr_mat, &mut curr_i, &mut curr_j, &mut aln)?;
} else {
self.retrace_mat1_mat3(&mut curr_mat, &mut curr_i, &mut curr_j, &mut aln)?;
}
}
let actual_aln = aln.reverse();
if has_ts {
Ok(TwitcherAlignmentCase::FoundTS {
alignment_with_ts: actual_aln,
cost_with_ts: best_result_mat3,
cost_without_ts: Some(best_result_mat1),
})
} else {
Ok(TwitcherAlignmentCase::NoTS {
alignment_without_ts: actual_aln,
cost_without_ts: best_result_mat1,
})
}
}
fn retrace_mat1_mat3(
&mut self,
curr_mat: &mut CurrentMatrix,
curr_i: &mut usize,
curr_j: &mut usize,
aln: &mut Alignment<AlignmentType>,
) -> anyhow::Result<()> {
use CurrentMatrix::{One, Two};
let curr_mat_ptr = if *curr_mat == One {
&self.mat1
} else {
&self.mat3
};
let cell = &curr_mat_ptr[(*curr_i, *curr_j)];
match cell.ty {
Some(Align::Match) => {
aln.push(AlignmentType::PrimaryMatch);
*curr_i -= 1;
*curr_j -= 1;
}
Some(Align::Substitution) => {
aln.push(AlignmentType::PrimarySubstitution);
*curr_i -= 1;
*curr_j -= 1;
}
Some(Align::Insertion) => {
aln.push(AlignmentType::PrimaryInsertion);
*curr_j -= 1;
}
Some(Align::Deletion) => {
aln.push(AlignmentType::PrimaryDeletion);
*curr_i -= 1;
}
Some(j @ (Align::MatchJump(new_i, new_j) | Align::SubstitutionJump(new_i, new_j))) => {
if matches!(j, Align::MatchJump(..)) {
aln.push(AlignmentType::PrimaryMatch);
} else {
aln.push(AlignmentType::PrimarySubstitution);
}
let len = self.mat2[(new_i, new_j)].ts_len;
let entrance_i = new_i + len as usize - 1;
let entrance_j = new_j - len as usize + 1;
let entrance_jump = &self.mat2[(entrance_i, entrance_j)].cell;
let Some(
Align::MatchJump(mat1_entrance_i, ..)
| Align::SubstitutionJump(mat1_entrance_i, ..),
) = entrance_jump.ty
else {
anyhow::bail!(
"We should have found the entrance jump in matrix 2 at i={entrance_i}, j={entrance_j}",
)
};
let ref_diff = isize::try_from(*curr_i)? - isize::try_from(mat1_entrance_i)? - 1;
aln.push(AlignmentType::TemplateSwitchExit {
anti_primary_gap: ref_diff,
});
*curr_mat = Two;
*curr_i = new_i;
*curr_j = new_j;
}
None => anyhow::bail!("Implementation error"),
}
Ok(())
}
fn retrace_mat2(
&mut self,
curr_mat: &mut CurrentMatrix,
curr_i: &mut usize,
curr_j: &mut usize,
aln: &mut Alignment<AlignmentType>,
) -> anyhow::Result<()> {
use CurrentMatrix::One;
let rev_cell = &self.mat2[(*curr_i, *curr_j)];
if aln.inner_mut().is_empty() {
// We are starting from matrix Two, i.e. we need to add the exit that is at the end of the alignment.
let len = rev_cell.ts_len;
let entrance_i = *curr_i + len as usize - 1;
let entrance_j = *curr_j - len as usize + 1;
let entrance_jump = &self.mat2[(entrance_i, entrance_j)].cell;
let Some(
Align::MatchJump(mat1_entrance_i, ..)
| Align::SubstitutionJump(mat1_entrance_i, ..),
) = entrance_jump.ty
else {
anyhow::bail!(
"We should have found the entrance jump in matrix 2 at i={entrance_i}, j={entrance_j}",
)
};
let ref_diff = isize::try_from(self.ranges.reference_range().len() + 1)?
- isize::try_from(mat1_entrance_i)?
- 1;
aln.push(AlignmentType::TemplateSwitchExit {
anti_primary_gap: ref_diff,
});
}
match rev_cell.cell.ty {
Some(Align::Match) => {
aln.push(AlignmentType::SecondaryMatch);
*curr_i += 1;
*curr_j -= 1;
}
Some(Align::Substitution) => {
aln.push(AlignmentType::SecondarySubstitution);
*curr_i += 1;
*curr_j -= 1;
}
Some(Align::Insertion) => {
aln.push(AlignmentType::SecondaryInsertion);
*curr_j -= 1;
}
Some(Align::Deletion) => {
aln.push(AlignmentType::SecondaryDeletion);
*curr_i += 1;
}
Some(j @ (Align::MatchJump(new_i, new_j) | Align::SubstitutionJump(new_i, new_j))) => {
if matches!(j, Align::MatchJump(..)) {
aln.push(AlignmentType::SecondaryMatch);
} else {
aln.push(AlignmentType::SecondarySubstitution);
}
aln.push(AlignmentType::TemplateSwitchEntrance { first_offset: isize::try_from(*curr_i)? - isize::try_from(new_i)?, equal_cost_range: EqualCostRange { min_start: 0, max_start: 0, min_end: 0, max_end: 0 }, primary: TemplateSwitchPrimary::Query, secondary: lib_tsalign::a_star_aligner::template_switch_distance::TemplateSwitchSecondary::Reference, direction: lib_tsalign::a_star_aligner::template_switch_distance::TemplateSwitchDirection::Reverse });
// ts_entrance_ref_ix
*curr_mat = One;
*curr_i = new_i;
*curr_j = new_j;
}
None => anyhow::bail!("Implementation error"),
}
Ok(())
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum CurrentMatrix {
One,
Two,
Three,
}
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()))
}
#[derive(Debug, PartialEq, Eq, Clone, Copy)]
enum Align {
Match,
Substitution,
Deletion,
Insertion,
MatchJump(usize, usize),
SubstitutionJump(usize, usize),
}
#[derive(Debug, PartialEq, Eq)]
struct Cell {
score: U64Cost,
ty: Option<Align>,
}
impl Display for Cell {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let op = match self.ty {
Some(Align::Match) => "Mat",
Some(Align::Substitution) => "Sub",
Some(Align::Deletion) => "Del",
Some(Align::Insertion) => "Ins",
Some(Align::MatchJump(_, _)) => "MaJ",
Some(Align::SubstitutionJump(_, _)) => "SuJ",
None => "???",
};
write!(f, "{:>3} ({op})", self.score.as_u64())
}
}
impl Default for Cell {
fn default() -> Self {
Self {
score: U64Cost::from_usize(0),
ty: None,
}
}
}
#[derive(Debug, Default, PartialEq, Eq)]
struct RevCell {
cell: Cell,
ts_len: u16,
}
impl Display for RevCell {
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 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
}
}
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::*;
#[test]
fn fpa_basic() {
let r = b"ACGTA";
let q = b"AACGTTA";
let fpa = FourPointAligner::default();
let r = fpa.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_complete(r.len(), q.len()),
);
r.unwrap();
}
#[test]
fn fpa_with_ranges() {
let r = b"TTTTTACGTATTTTT";
let q = b"TTTTTACGAAAAAAATTTTT";
let fpa = FourPointAligner::default();
let r = fpa.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(5, 5),
AlignmentCoordinates::new(10, 15),
),
);
r.unwrap();
}
#[test]
fn fpa_with_cursed_ranges() {
let r = b"TTTTTACGTATTTTT";
let q = b"TTTTTACGAAAAAAATTTTT";
let fpa = FourPointAligner::default();
let r = fpa.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(0, 18),
AlignmentCoordinates::new(1, 20),
),
);
r.unwrap();
}
#[test]
fn fpa_vcf() {
let r = b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA";
let q = b"AAAAAAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAAA";
let fpa = FourPointAligner::default();
let r = fpa.align(
Arc::from(&r[..]),
Arc::from(&q[..]),
AlignmentRange::new_offset_limit(
AlignmentCoordinates::new(15, 15),
AlignmentCoordinates::new(26, 25),
),
);
dbg!(&r);
r.unwrap();
}
#[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!(res.has_ts());
}
}