26 #include <range/v3/view/zip.hpp> 31 struct view_equality_fn
34 template <std::ranges::ForwardRange rng1_type, std::ranges::ForwardRange rng2_type>
35 constexpr
bool operator()(rng1_type && rng1, rng2_type && rng2)
const 84 template <
typename reference_
char_type,
typename query_
char_type>
86 requires std::detail::WeaklyEqualityComparableWith<reference_char_type, gap> &&
87 std::detail::WeaklyEqualityComparableWith<query_char_type, gap>
89 char compare_aligned_values(reference_char_type
const reference_char,
90 query_char_type
const query_char,
91 bool const extended_cigar)
93 return (reference_char == gap{})
94 ? (query_char == gap{})
97 : (query_char == gap{})
100 ? (query_char == reference_char)
137 template <std::ranges::ForwardRange ref_seq_type, std::ranges::ForwardRange query_seq_type>
139 requires std::detail::WeaklyEqualityComparableWith<gap, reference_t<ref_seq_type>> &&
140 std::detail::WeaklyEqualityComparableWith<gap, reference_t<query_seq_type>>
142 std::string get_cigar_string(ref_seq_type && ref_seq,
143 query_seq_type && query_seq,
144 uint32_t
const query_start_pos = 0,
145 uint32_t
const query_end_pos = 0,
146 bool const extended_cigar =
false)
148 if (ref_seq.size() != query_seq.size())
158 result << query_start_pos <<
'S';
163 char tmp_char{compare_aligned_values(ref_seq[0], query_seq[0], extended_cigar)};
164 size_t tmp_length{0};
169 char next_op = compare_aligned_values(std::get<0>(column), std::get<1>(column), extended_cigar);
171 if (tmp_char == next_op)
177 result << tmp_length << tmp_char;
183 result << tmp_length << tmp_char;
187 result << query_end_pos <<
'S';
231 template <TupleLike alignment_type>
233 requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2
235 std::string get_cigar_string(alignment_type && alignment,
236 uint32_t
const query_start_pos = 0,
237 uint32_t
const query_end_pos = 0,
238 bool const extended_cigar =
false)
240 return get_cigar_string(get<0>(alignment), get<1>(alignment),
241 query_start_pos, query_end_pos, extended_cigar);
282 template <TupleLike alignment_type>
284 requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2
287 uint32_t
const query_start_pos = 0,
288 uint32_t
const query_end_pos = 0,
289 bool const extended_cigar =
false)
291 auto & ref_seq = get<0>(alignment);
292 auto & query_seq = get<1>(alignment);
294 if (ref_seq.size() != query_seq.size())
304 result.emplace_back(
'S', query_start_pos);
309 char tmp_char{compare_aligned_values(ref_seq[0], query_seq[0], extended_cigar)};
310 size_t tmp_length{0};
315 char next_op = compare_aligned_values(std::get<0>(column), std::get<1>(column), extended_cigar);
317 if (tmp_char == next_op)
323 result.emplace_back(tmp_char, tmp_length);
330 result.emplace_back(tmp_char, tmp_length);
334 result.emplace_back(
'S', query_end_pos);
353 template <std::ranges::InputRange cigar_input_type>
355 parse_cigar(cigar_input_type && cigar_input)
358 size_t sc_begin_count{};
359 size_t sc_end_count{};
362 size_t cigar_count{}, ref_length{}, seq_length{};
364 auto update_lengths_fn = [&ref_length, &seq_length, &cigar_op, &cigar_count] ()
366 if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
368 ref_length += cigar_count;
369 seq_length += cigar_count;
371 else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op))
373 ref_length += cigar_count;
375 else if (is_char<'I'>(cigar_op))
377 seq_length += cigar_count;
381 if (is_char<'P'>(cigar_op))
382 throw format_error{
"We do currently not support cigar operation 'P'."};
399 if (is_char<'H'>(cigar_op))
403 buffer_end = buffer_end2;
411 throw format_error{
"Corrupted cigar string encountered"};
413 if (is_char<'S'>(cigar_op))
415 sc_begin_count = cigar_count;
420 operations.push_back({cigar_op, cigar_count});
432 throw format_error{
"Corrupted cigar string encountered"};
434 if (is_char<'S'>(cigar_op))
436 sc_end_count = cigar_count;
437 return {operations, ref_length, seq_length, sc_begin_count, sc_end_count};
440 operations.push_back({cigar_op, cigar_count});
442 return {operations, ref_length, seq_length, sc_begin_count, sc_end_count};
460 template <std::ranges::InputRange cigar_input_type>
461 auto parse_binary_cigar(cigar_input_type && cigar_input, uint16_t n_cigar_op)
464 size_t sc_begin_count{};
465 size_t sc_end_count{};
467 size_t cigar_count{},ref_length{}, seq_length{};
469 static constexpr
char const * CIGAR_MAPPING =
"MIDNSHP=X*******";
470 static constexpr uint32_t CIGAR_MASK = 0x0f;
472 auto update_lengths_fn = [&ref_length, &seq_length, &cigar_op, &cigar_count] ()
474 if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
476 ref_length += cigar_count;
477 seq_length += cigar_count;
479 else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op))
481 ref_length += cigar_count;
483 else if (is_char<'I'>(cigar_op))
485 seq_length += cigar_count;
489 if (is_char<'P'>(cigar_op))
490 throw format_error{
"We do currently not support cigar operation 'P'."};
497 return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
500 assert((o_and_c & CIGAR_MASK) <= 8u);
501 cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
502 cigar_count = o_and_c >> 4;
505 if (is_char<'H'>(cigar_op))
508 return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
511 assert((o_and_c & CIGAR_MASK) <= 8u);
512 cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
513 cigar_count = o_and_c >> 4;
517 if (is_char<'S'>(cigar_op))
519 sc_begin_count = cigar_count;
524 operations.push_back({cigar_op, cigar_count});
529 while (n_cigar_op > 0)
532 assert((o_and_c & CIGAR_MASK) <= 8u);
533 cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
534 cigar_count = o_and_c >> 4;
536 if (is_char<'S'>(cigar_op))
538 sc_end_count = cigar_count;
539 return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
543 operations.push_back({cigar_op, cigar_count});
546 return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
572 template <TupleLike alignment_type>
574 requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2 &&
575 detail::all_satisfy_aligned_seq<detail::tuple_type_list_t<alignment_type>>
583 for (
auto [cigar_op, cigar_count] : cigar)
585 if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
590 else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op))
593 current_read_pos =
insert_gap(get<1>(alignment), current_read_pos, cigar_count);
597 else if (is_char<'I'>(cigar_op))
600 current_ref_pos =
insert_gap(get<0>(alignment), current_ref_pos, cigar_count);
606 if (is_char<'P'>(cigar_op))
607 throw format_error{
"We do currently not support cigar operation 'P'."};
615 struct access_restrictor_fn
618 template <
typename chr_t>
619 [[noreturn]] chr_t operator()(chr_t)
const 621 throw std::logic_error{
"Access is not allowed because there is no sequence information."};
::ranges::equal equal
Alias for ranges::equal. Determines if two sets of elements are the same.
Definition: algorithm:54
::ranges::distance distance
Alias for ranges::distance. Returns the number of hops from first to last.
Definition: iterator:321
::ranges::next next
Alias for ranges::next. Returns the nth successor of the given iterator.
Definition: iterator:331
aligned_seq_t::iterator insert_gap(aligned_seq_t &aligned_seq, typename aligned_seq_t::const_iterator pos_it)
An implementation of seqan3::AlignedSequence::insert_gap for sequence containers. ...
Definition: aligned_sequence_concept.hpp:229
Includes the AlignedSequence and the related insert_gap and erase_gap functions to enable stl contain...
constexpr auto zip
A range adaptor that transforms a tuple of range into a range of tuples.
Definition: ranges:948
::ranges::copy_n copy_n
Alias for ranges::copy_n. Copies a range of exactly n elements to a new location. ...
Definition: algorithm:49
Provides std::from_chars and std::to_chars if not defined in the stl <charconv> header.
Provides seqan3::view::take_until and seqan3::view::take_until_or_throw.
Provides seqan3::TupleLike.
auto constexpr is_digit
Checks whether c is a digital character.
Definition: predicate.hpp:287
Adaptations of concepts from the Ranges TS.
::ranges::begin begin
Alias for ranges::begin. Returns an iterator to the beginning of a range.
Definition: ranges:174
::ranges::advance advance
Alias for ranges::advance. Advances the iterator by the given distance.
Definition: iterator:316
::ranges::copy copy
Alias for ranges::copy. Copies a range of elements to a new location.
Definition: algorithm:44
constexpr auto single_pass_input
A view adapter that decays most of the range properties and adds single pass behavior.
Definition: single_pass_input.hpp:381
Definition: aligned_sequence_concept.hpp:35
auto const get
A view calling std::get on each element in a range.
Definition: get.hpp:66
auto constexpr take_until_or_throw
A view adaptor that returns elements from the underlying range until the functor evaluates to true (t...
Definition: take_until.hpp:613
Adaptations of algorithms from the Ranges TS.
::ranges::end end
Alias for ranges::end. Returns an iterator to the end of a range.
Definition: ranges:179