83template <
typename reference_type,
typename sequence_type>
85 reference_type
const & reference,
86 uint32_t
const zero_based_reference_start_position,
87 sequence_type
const & query)
89 if (cigar_vector.
empty())
90 throw std::logic_error{
"An empty CIGAR is not a valid alignment representation."};
95 uint32_t reference_length{0};
96 uint32_t query_length{0};
98 for (
auto [cigar_count, cigar_operation] : cigar_vector)
100 if ((
'M'_cigar_operation == cigar_operation) || (
'='_cigar_operation == cigar_operation)
101 || (
'X'_cigar_operation == cigar_operation))
103 reference_length += cigar_count;
104 query_length += cigar_count;
106 else if (
'D'_cigar_operation == cigar_operation)
108 reference_length += cigar_count;
110 else if (
'I'_cigar_operation == cigar_operation)
112 query_length += cigar_count;
116 if (
static_cast<size_t>(zero_based_reference_start_position + reference_length) >
std::ranges::size(reference))
117 throw std::logic_error{
"The CIGAR string indicates a reference length of at least "
118 +
std::to_string(zero_based_reference_start_position + reference_length)
119 +
", but the supplied reference sequence is only of size"
124 uint32_t soft_clipping_start{0};
125 uint32_t soft_clipping_end{0};
128 auto soft_clipping_at = [&cigar_vector](
size_t const index)
130 return cigar_vector[index] ==
'S'_cigar_operation;
133 auto hard_clipping_at = [&](
size_t const index)
135 return cigar_vector[index] ==
'H'_cigar_operation;
138 auto vector_size_at_least = [&](
size_t const min_size)
140 return cigar_vector.
size() >= min_size;
143 auto cigar_count_at = [&](
size_t const index)
145 return get<0>(cigar_vector[index]);
150 if (soft_clipping_at(0))
151 soft_clipping_start = cigar_count_at(0);
152 else if (vector_size_at_least(2) && hard_clipping_at(0) && soft_clipping_at(1))
153 soft_clipping_start = cigar_count_at(1);
157 auto const last_index = cigar_vector.
size() - 1;
158 auto const second_last_index = last_index - 1;
160 if (vector_size_at_least(2) && soft_clipping_at(last_index))
161 soft_clipping_end = cigar_count_at(last_index);
162 else if (vector_size_at_least(3) && hard_clipping_at(last_index) && soft_clipping_at(second_last_index))
163 soft_clipping_end = cigar_count_at(second_last_index);
165 if (soft_clipping_start + query_length + soft_clipping_end !=
std::ranges::size(query))
166 throw std::logic_error{
"The CIGAR string indicates a query/read sequence length of "
167 +
std::to_string(soft_clipping_start + query_length + soft_clipping_end)
168 +
", but the supplied query/read sequence is of size"
182 zero_based_reference_start_position + reference_length));
184 assign_unaligned(get<1>(
alignment), query |
views::slice(soft_clipping_start, soft_clipping_start + query_length));
192 for (
auto [cigar_count, cigar_operation] : cigar_vector)
194 if ((
'M'_cigar_operation == cigar_operation) || (
'='_cigar_operation == cigar_operation)
195 || (
'X'_cigar_operation == cigar_operation))
197 std::ranges::advance(current_ref_pos, cigar_count);
198 std::ranges::advance(current_read_pos, cigar_count);
200 else if ((
'D'_cigar_operation == cigar_operation) || (
'N'_cigar_operation == cigar_operation))
203 current_read_pos = get<1>(
alignment).insert_gap(current_read_pos, cigar_count);
205 std::ranges::advance(current_ref_pos, cigar_count);
207 else if ((
'I'_cigar_operation == cigar_operation))
209 current_ref_pos = get<0>(
alignment).insert_gap(current_ref_pos, cigar_count);
211 std::ranges::advance(current_read_pos, cigar_count);
213 else if ((
'P'_cigar_operation == cigar_operation))
215 current_ref_pos = get<0>(
alignment).insert_gap(current_ref_pos, cigar_count);
218 current_read_pos = get<1>(
alignment).insert_gap(current_read_pos, cigar_count);
230template <
typename reference_type,
typename sequence_type>
232 reference_type
const & reference,
233 uint32_t
const zero_based_reference_start_position,
234 sequence_type
const & query)
Includes the aligned_sequence and the related insert_gap and erase_gap functions to enable stl contai...
Provides the seqan3::cigar alphabet.
A gap decorator allows the annotation of sequences with gap symbols while leaving the underlying sequ...
Definition: gap_decorator.hpp:81
Provides seqan3::gap_decorator.
auto alignment_from_cigar(std::vector< cigar > const &cigar_vector, reference_type const &reference, uint32_t const zero_based_reference_start_position, sequence_type const &query)
Construct an alignment from a CIGAR string and the corresponding sequences.
Definition: alignment_from_cigar.hpp:84
constexpr std::tuple< std::vector< cigar >, int32_t, int32_t > parse_cigar(cigar_input_type &&cigar_input)
Parses a cigar string into a vector of operation-count pairs (e.g. (M, 3)).
Definition: io/sam_file/detail/cigar.hpp:94
@ alignment
The (pairwise) alignment stored in an object that models seqan3::detail::pairwise_alignment.
constexpr size_t size
The size of a type pack.
Definition: type_pack/traits.hpp:146
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition: slice.hpp:178
Auxiliary functions for the SAM IO.
The main SeqAn3 namespace.
Definition: aligned_sequence_concept.hpp:29
constexpr auto const & get(configuration< configs_t... > const &config) noexcept
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition: configuration.hpp:415
Provides seqan3::views::slice.