use crate::{
alignment::{AlignmentIter, AlignmentStates, pairwise_align_with},
data::{cigar::Ciglet, views::IndexAdjustable},
};
use std::ops::Range;
#[derive(Clone, Eq, PartialEq, Debug)]
pub enum MaybeAligned<T> {
Some(T),
Overflowed,
Unmapped,
}
impl<T> MaybeAligned<T> {
#[inline]
#[must_use]
pub fn unwrap(self) -> T {
match self {
MaybeAligned::Some(aln) => aln,
MaybeAligned::Overflowed => panic!("Alignment score overflowed!"),
MaybeAligned::Unmapped => panic!("Sequence could not be mapped!"),
}
}
#[inline]
#[must_use]
pub fn get(self) -> Option<T> {
match self {
MaybeAligned::Some(aln) => Some(aln),
_ => None,
}
}
#[inline]
#[must_use]
pub const fn as_ref(&self) -> MaybeAligned<&T> {
match *self {
MaybeAligned::Some(ref x) => MaybeAligned::Some(x),
MaybeAligned::Overflowed => MaybeAligned::Overflowed,
MaybeAligned::Unmapped => MaybeAligned::Unmapped,
}
}
#[inline]
#[must_use]
pub const fn as_mut(&mut self) -> MaybeAligned<&mut T> {
match *self {
MaybeAligned::Some(ref mut x) => MaybeAligned::Some(x),
MaybeAligned::Overflowed => MaybeAligned::Overflowed,
MaybeAligned::Unmapped => MaybeAligned::Unmapped,
}
}
#[inline]
#[must_use]
pub fn or_else_overflowed(self, f: impl FnOnce() -> Self) -> Self {
if let MaybeAligned::Overflowed = self { f() } else { self }
}
#[inline]
pub fn map<U>(self, f: impl FnOnce(T) -> U) -> MaybeAligned<U> {
match self {
MaybeAligned::Some(value) => MaybeAligned::Some(f(value)),
MaybeAligned::Overflowed => MaybeAligned::Overflowed,
MaybeAligned::Unmapped => MaybeAligned::Unmapped,
}
}
#[inline]
#[must_use]
pub fn and_then<U, F>(self, f: F) -> MaybeAligned<U>
where
F: FnOnce(T) -> MaybeAligned<U>, {
match self {
MaybeAligned::Some(x) => f(x),
MaybeAligned::Overflowed => MaybeAligned::Overflowed,
MaybeAligned::Unmapped => MaybeAligned::Unmapped,
}
}
}
impl<T> MaybeAligned<Alignment<T>>
where
T: Copy,
{
#[inline]
#[must_use]
pub fn to_reverse(&self) -> Self {
self.as_ref().map(Alignment::to_reverse)
}
#[inline]
pub fn make_reverse(&mut self) {
self.as_mut().map(Alignment::make_reverse);
}
}
pub(crate) struct ScoreIndices<T> {
pub score: T,
pub ref_idx: usize,
pub query_idx: usize,
}
impl<T> ScoreIndices<T> {
#[inline]
#[must_use]
pub fn into_score_starts(self) -> ScoreStarts<T> {
ScoreStarts {
score: self.score,
ref_start: self.ref_idx,
query_start: self.query_idx,
}
}
#[inline]
#[must_use]
pub fn into_score_ends(self) -> ScoreEnds<T> {
ScoreEnds {
score: self.score,
ref_end: self.ref_idx,
query_end: self.query_idx,
}
}
}
pub struct ScoreStarts<T> {
pub score: T,
pub ref_start: usize,
pub query_start: usize,
}
pub struct ScoreEnds<T> {
pub score: T,
pub ref_end: usize,
pub query_end: usize,
}
pub struct ScoreAndRanges<T> {
pub score: T,
pub ref_range: Range<usize>,
pub query_range: Range<usize>,
}
#[non_exhaustive]
#[derive(Clone, Eq, PartialEq, Debug, Default)]
pub struct Alignment<T> {
pub score: T,
pub ref_range: Range<usize>,
pub query_range: Range<usize>,
pub states: AlignmentStates,
pub ref_len: usize,
pub query_len: usize,
}
impl<T> Alignment<T> {
#[inline]
#[must_use]
pub fn new_global(score: T, states: AlignmentStates, ref_len: usize, query_len: usize) -> Self {
Self {
score,
ref_range: 0..ref_len,
query_range: 0..query_len,
states,
ref_len,
query_len,
}
}
#[inline]
#[must_use]
pub fn get_aligned_seqs(&self, reference: &[u8], query: &[u8]) -> (Vec<u8>, Vec<u8>) {
pairwise_align_with(reference, query, &self.states, self.ref_range.start)
}
#[inline]
#[must_use]
pub fn get_aligned_iter<'a>(
&'a self, reference: &'a [u8], query: &'a [u8],
) -> AlignmentIter<'a, std::iter::Copied<std::slice::Iter<'a, Ciglet>>> {
AlignmentIter::new(reference, query, &self.states, self.ref_range.start)
}
pub fn get_aligned_query(&self, query: &[u8]) -> Vec<u8> {
let mut query_index = 0;
let mut query_aln = Vec::with_capacity(query.len() + (self.ref_range.len() / 2));
for Ciglet { inc, op } in &self.states {
match op {
b'M' | b'=' | b'X' | b'I' => {
query_aln.extend_from_slice(&query[query_index..query_index + inc]);
query_index += inc;
}
b'D' => {
query_aln.extend(std::iter::repeat_n(b'-', inc));
}
b'S' => query_index += inc,
b'N' => {
query_aln.extend(std::iter::repeat_n(b'N', inc));
}
b'H' | b'P' => {}
_ => panic!("CIGAR op '{op}' not supported.\n"),
}
}
query_aln
}
}
impl<T: Copy> Alignment<T> {
#[must_use]
pub fn invert(&self) -> Self {
let mut states = AlignmentStates::new();
states.soft_clip(self.ref_range.start);
let inverted_ciglets = self.states.0.iter().filter_map(|&ciglet| match ciglet.op {
b'S' | b'H' => None,
b'D' => Some(Ciglet {
inc: ciglet.inc,
op: b'I',
}),
b'I' => Some(Ciglet {
inc: ciglet.inc,
op: b'D',
}),
_ => Some(ciglet),
});
states.extend_from_ciglets(inverted_ciglets);
states.soft_clip(self.ref_len - self.ref_range.end);
Self {
score: self.score,
ref_range: self.query_range.clone(),
query_range: self.ref_range.clone(),
states,
ref_len: self.query_len,
query_len: self.ref_len,
}
}
#[must_use]
pub fn to_reverse(&self) -> Self {
let ref_range = (self.ref_len - self.ref_range.end)..(self.ref_len - self.ref_range.start);
let query_range = (self.query_len - self.query_range.end)..(self.query_len - self.query_range.start);
let states = self.states.to_reverse();
Alignment {
score: self.score,
ref_range,
query_range,
states,
ref_len: self.ref_len,
query_len: self.query_len,
}
}
pub fn make_reverse(&mut self) {
self.ref_range = (self.ref_len - self.ref_range.end)..(self.ref_len - self.ref_range.start);
self.query_range = (self.query_len - self.query_range.end)..(self.query_len - self.query_range.start);
self.states.make_reverse();
}
#[inline]
#[must_use]
pub fn to_verbose_sequence_matching(&self, reference: &[u8], query: &[u8]) -> Self {
Self {
score: self.score,
ref_range: self.ref_range.clone(),
query_range: self.query_range.clone(),
states: self
.states
.to_verbose_sequence_matching(reference, query, self.ref_range.start),
ref_len: self.ref_len,
query_len: self.query_len,
}
}
pub fn slice_to_ref_range(&self, ref_range: Range<usize>) -> Option<Self> {
let rel_ref_range = if ref_range.start >= self.ref_range.start {
ref_range.sub(self.ref_range.start)
} else {
return None;
};
let mut query_idx = 0;
let mut ref_idx = 0;
let mut ciglets = self.states.iter().copied();
let mut states = AlignmentStates::new();
let Some(mut first_ciglet) = ciglets.find_map(|ciglet| {
(query_idx, ref_idx) = (query_idx, ref_idx).increment_idxs_by(ciglet);
let num_in_range = ref_idx.saturating_sub(rel_ref_range.start);
(num_in_range > 0).then_some(Ciglet {
op: ciglet.op,
inc: num_in_range,
})
}) else {
if ref_range.is_empty() && (self.ref_range.start..=self.ref_range.end).contains(&ref_range.start) {
let mut states = AlignmentStates::with_capacity(1);
states.soft_clip(self.query_len);
return Some(Alignment {
score: self.score,
ref_range: ref_range.clone(),
query_range: 0..0,
states,
ref_len: self.ref_len,
query_len: self.query_len,
});
}
if !self.ref_range.is_empty() {
debug_assert!(ref_range.start >= self.ref_range.end);
}
return None;
};
let (query_start, ref_start) = (query_idx, ref_idx).decrement_idxs_by(first_ciglet);
debug_assert_eq!(ref_start, rel_ref_range.start);
states.soft_clip(query_start);
if ref_idx >= rel_ref_range.end {
first_ciglet.inc = first_ciglet.inc.saturating_sub(ref_idx - rel_ref_range.end);
states.add_ciglet(first_ciglet);
let (query_end, ref_end) = (query_start, ref_start).increment_idxs_by(first_ciglet);
states.soft_clip(self.query_len - query_end);
debug_assert_eq!(ref_end, rel_ref_range.end);
return Some(Self {
score: self.score,
ref_range,
query_range: query_start..query_end,
states,
ref_len: self.ref_len,
query_len: self.query_len,
});
}
states.add_ciglet(first_ciglet);
for mut ciglet in ciglets {
let (new_query_idx, new_ref_idx) = (query_idx, ref_idx).increment_idxs_by(ciglet);
if let Some(num_past) = new_ref_idx.checked_sub(rel_ref_range.end) {
ciglet.inc = ciglet.inc.saturating_sub(num_past);
states.add_ciglet(ciglet);
let (query_end, ref_end) = (query_idx, ref_idx).increment_idxs_by(ciglet);
states.soft_clip(self.query_len - query_end);
debug_assert!(ciglet.inc > 0);
debug_assert_eq!(ref_end, rel_ref_range.end);
return Some(Self {
score: self.score,
ref_range,
query_range: query_start..query_end,
states,
ref_len: self.ref_len,
query_len: self.query_len,
});
}
states.add_ciglet(ciglet);
(query_idx, ref_idx) = (new_query_idx, new_ref_idx);
}
if ref_idx < rel_ref_range.end {
debug_assert!(ref_range.end > self.ref_range.end);
return None;
}
debug_assert_eq!(ref_idx, rel_ref_range.end);
let query_range = query_start..query_idx;
states.soft_clip(self.query_len - query_idx);
Some(Self {
score: self.score,
ref_range,
query_range,
states,
ref_len: self.ref_len,
query_len: self.query_len,
})
}
}
pub(crate) trait AlignmentIndices {
fn increment_idxs_by(self, ciglet: Ciglet) -> Self;
fn decrement_idxs_by(self, ciglet: Ciglet) -> Self;
}
impl AlignmentIndices for (usize, usize) {
#[inline]
fn increment_idxs_by(self, ciglet: Ciglet) -> Self {
let (query_idx, ref_idx) = self;
match ciglet.op {
b'M' | b'=' | b'X' => (query_idx + ciglet.inc, ref_idx + ciglet.inc),
b'D' | b'N' => (query_idx, ref_idx + ciglet.inc),
b'I' | b'S' => (query_idx + ciglet.inc, ref_idx),
_ => (query_idx, ref_idx),
}
}
#[inline]
fn decrement_idxs_by(self, ciglet: Ciglet) -> Self {
let (query_idx, ref_idx) = self;
match ciglet.op {
b'M' | b'=' | b'X' => (query_idx - ciglet.inc, ref_idx - ciglet.inc),
b'D' | b'N' => (query_idx, ref_idx - ciglet.inc),
b'I' | b'S' => (query_idx - ciglet.inc, ref_idx),
_ => (query_idx, ref_idx),
}
}
}