82template <
typename reference_type,
typename sequence_type>
84 reference_type
const & reference,
85 uint32_t
const zero_based_reference_start_position,
86 sequence_type
const & query)
88 if (cigar_vector.
empty())
89 throw std::logic_error{
"An empty CIGAR is not a valid alignment representation."};
94 uint32_t reference_length{0};
95 uint32_t query_length{0};
97 for (
auto [cigar_count, cigar_operation] : cigar_vector)
99 if ((
'M'_cigar_operation == cigar_operation) || (
'='_cigar_operation == cigar_operation)
100 || (
'X'_cigar_operation == cigar_operation))
102 reference_length += cigar_count;
103 query_length += cigar_count;
105 else if (
'D'_cigar_operation == cigar_operation)
107 reference_length += cigar_count;
109 else if (
'I'_cigar_operation == cigar_operation)
111 query_length += cigar_count;
115 if (
static_cast<size_t>(zero_based_reference_start_position + reference_length) >
std::ranges::size(reference))
116 throw std::logic_error{
"The CIGAR string indicates a reference length of at least "
117 +
std::to_string(zero_based_reference_start_position + reference_length)
118 +
", but the supplied reference sequence is only of size"
123 uint32_t soft_clipping_start{0};
124 uint32_t soft_clipping_end{0};
127 auto soft_clipping_at = [&cigar_vector](
size_t const index)
129 return cigar_vector[index] ==
'S'_cigar_operation;
132 auto hard_clipping_at = [&](
size_t const index)
134 return cigar_vector[index] ==
'H'_cigar_operation;
137 auto vector_size_at_least = [&](
size_t const min_size)
139 return cigar_vector.
size() >= min_size;
142 auto cigar_count_at = [&](
size_t const index)
144 return get<0>(cigar_vector[index]);
149 if (soft_clipping_at(0))
150 soft_clipping_start = cigar_count_at(0);
151 else if (vector_size_at_least(2) && hard_clipping_at(0) && soft_clipping_at(1))
152 soft_clipping_start = cigar_count_at(1);
156 auto const last_index = cigar_vector.
size() - 1;
157 auto const second_last_index = last_index - 1;
159 if (vector_size_at_least(2) && soft_clipping_at(last_index))
160 soft_clipping_end = cigar_count_at(last_index);
161 else if (vector_size_at_least(3) && hard_clipping_at(last_index) && soft_clipping_at(second_last_index))
162 soft_clipping_end = cigar_count_at(second_last_index);
164 if (soft_clipping_start + query_length + soft_clipping_end !=
std::ranges::size(query))
165 throw std::logic_error{
"The CIGAR string indicates a query/read sequence length of "
166 +
std::to_string(soft_clipping_start + query_length + soft_clipping_end)
167 +
", but the supplied query/read sequence is of size"
181 zero_based_reference_start_position + reference_length));
183 assign_unaligned(get<1>(
alignment), query |
views::slice(soft_clipping_start, soft_clipping_start + query_length));
191 for (
auto [cigar_count, cigar_operation] : cigar_vector)
193 if ((
'M'_cigar_operation == cigar_operation) || (
'='_cigar_operation == cigar_operation)
194 || (
'X'_cigar_operation == cigar_operation))
196 std::ranges::advance(current_ref_pos, cigar_count);
197 std::ranges::advance(current_read_pos, cigar_count);
199 else if ((
'D'_cigar_operation == cigar_operation) || (
'N'_cigar_operation == cigar_operation))
202 current_read_pos = get<1>(
alignment).insert_gap(current_read_pos, cigar_count);
204 std::ranges::advance(current_ref_pos, cigar_count);
206 else if ((
'I'_cigar_operation == cigar_operation))
208 current_ref_pos = get<0>(
alignment).insert_gap(current_ref_pos, cigar_count);
210 std::ranges::advance(current_read_pos, cigar_count);
212 else if ((
'P'_cigar_operation == cigar_operation))
214 current_ref_pos = get<0>(
alignment).insert_gap(current_ref_pos, cigar_count);
217 current_read_pos = get<1>(
alignment).insert_gap(current_read_pos, cigar_count);
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:83
@ 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
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.