Loading...
Searching...
No Matches
alignment_from_cigar.hpp
Go to the documentation of this file.
// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
#pragma once
#include <vector>
#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/alignment/decorator/gap_decorator.hpp>
#include <seqan3/alphabet/cigar/cigar.hpp>
#include <seqan3/io/sam_file/detail/cigar.hpp>
#include <seqan3/utility/views/slice.hpp>
namespace seqan3
{
template <typename reference_type, typename sequence_type>
inline 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)
{
if (cigar_vector.empty())
throw std::logic_error{"An empty CIGAR is not a valid alignment representation."};
// compute the length of the aligned region in the reference sequence
// -------------------------------------------------------------------------
// this requires a first stream over the cigar vector.
uint32_t reference_length{0};
uint32_t query_length{0};
for (auto [cigar_count, cigar_operation] : cigar_vector)
{
if (('M'_cigar_operation == cigar_operation) || ('='_cigar_operation == cigar_operation)
|| ('X'_cigar_operation == cigar_operation))
{
reference_length += cigar_count;
query_length += cigar_count;
}
else if ('D'_cigar_operation == cigar_operation)
{
reference_length += cigar_count;
}
else if ('I'_cigar_operation == cigar_operation)
{
query_length += cigar_count;
}
}
if (static_cast<size_t>(zero_based_reference_start_position + reference_length) > std::ranges::size(reference))
throw std::logic_error{"The CIGAR string indicates a reference length of at least "
+ std::to_string(zero_based_reference_start_position + reference_length)
+ ", but the supplied reference sequence is only of size"
+ std::to_string(std::ranges::size(reference)) + "."};
// Get soft clipping at the start and end of the CIGAR string
// -------------------------------------------------------------------------
uint32_t soft_clipping_start{0};
uint32_t soft_clipping_end{0};
// Checks whether the given index in the cigar vector is a soft clip.
auto soft_clipping_at = [&cigar_vector](size_t const index)
{
return cigar_vector[index] == 'S'_cigar_operation;
};
// Checks whether the given index in the cigar vector is a hard clip.
auto hard_clipping_at = [&](size_t const index)
{
return cigar_vector[index] == 'H'_cigar_operation;
};
// Checks whether the given cigar vector has at least min_size many elements.
auto vector_size_at_least = [&](size_t const min_size)
{
return cigar_vector.size() >= min_size;
};
// Returns the cigar count of the ith cigar element in the given cigar vector.
auto cigar_count_at = [&](size_t const index)
{
return get<0>(cigar_vector[index]);
};
// check for soft clipping at the first two positions
// cigar is non-empty, checked at the very beginning.
if (soft_clipping_at(0))
soft_clipping_start = cigar_count_at(0);
else if (vector_size_at_least(2) && hard_clipping_at(0) && soft_clipping_at(1))
soft_clipping_start = cigar_count_at(1);
// Check for soft clipping at the last two positions to validate the CIGAR string.
// Even if the two following arithmetics overflow, they are protected by the corresponding if expressions below.
auto const last_index = cigar_vector.size() - 1;
auto const second_last_index = last_index - 1;
if (vector_size_at_least(2) && soft_clipping_at(last_index))
soft_clipping_end = cigar_count_at(last_index);
else if (vector_size_at_least(3) && hard_clipping_at(last_index) && soft_clipping_at(second_last_index))
soft_clipping_end = cigar_count_at(second_last_index);
if (soft_clipping_start + query_length + soft_clipping_end != std::ranges::size(query))
throw std::logic_error{"The CIGAR string indicates a query/read sequence length of "
+ std::to_string(soft_clipping_start + query_length + soft_clipping_end)
+ ", but the supplied query/read sequence is of size"
+ std::to_string(std::ranges::size(query)) + "."};
// Assign the sequence to the alignment (a tuple of 2 gap decorators)
// -------------------------------------------------------------------------
using gapped_reference_type = gap_decorator<decltype(reference | views::slice(0, 0))>;
using gapped_sequence_type = gap_decorator<decltype(query | views::slice(0, 0))>;
using alignment_type = std::tuple<gapped_reference_type, gapped_sequence_type>;
alignment_type alignment{};
assign_unaligned(get<0>(alignment),
reference
| views::slice(zero_based_reference_start_position,
zero_based_reference_start_position + reference_length));
// query_length already accounts for soft clipping at begin and end
assign_unaligned(get<1>(alignment), query | views::slice(soft_clipping_start, soft_clipping_start + query_length));
// Insert gaps into the alignment based on the cigar vector
// -------------------------------------------------------------------------
using std::get;
auto current_ref_pos = std::ranges::begin(get<0>(alignment));
auto current_read_pos = std::ranges::begin(get<1>(alignment));
for (auto [cigar_count, cigar_operation] : cigar_vector)
{
if (('M'_cigar_operation == cigar_operation) || ('='_cigar_operation == cigar_operation)
|| ('X'_cigar_operation == cigar_operation))
{
std::ranges::advance(current_ref_pos, cigar_count);
std::ranges::advance(current_read_pos, cigar_count);
}
else if (('D'_cigar_operation == cigar_operation) || ('N'_cigar_operation == cigar_operation))
{
// insert gaps into query
current_read_pos = get<1>(alignment).insert_gap(current_read_pos, cigar_count);
++current_read_pos;
std::ranges::advance(current_ref_pos, cigar_count);
}
else if (('I'_cigar_operation == cigar_operation)) // Insert gaps into ref
{
current_ref_pos = get<0>(alignment).insert_gap(current_ref_pos, cigar_count);
++current_ref_pos;
std::ranges::advance(current_read_pos, cigar_count);
}
else if (('P'_cigar_operation == cigar_operation)) // skip padding
{
current_ref_pos = get<0>(alignment).insert_gap(current_ref_pos, cigar_count);
++current_ref_pos;
current_read_pos = get<1>(alignment).insert_gap(current_read_pos, cigar_count);
++current_read_pos;
}
// S and H are ignored as they are handled by cropping the sequence
}
return alignment;
}
template <typename reference_type, typename sequence_type>
inline auto alignment_from_cigar(std::string const & cigar_string,
reference_type const & reference,
uint32_t const zero_based_reference_start_position,
sequence_type const & query)
{
alignment_from_cigar(detail::parse_cigar(cigar_string), reference, zero_based_reference_start_position, query);
}
} // namespace seqan3
aligned_sequence_concept.hpp
Includes the aligned_sequence and the related insert_gap and erase_gap functions to enable stl contai...
cigar.hpp
Provides the seqan3::cigar alphabet.
std::string
seqan3::gap_decorator
A gap decorator allows the annotation of sequences with gap symbols while leaving the underlying sequ...
Definition gap_decorator.hpp:78
std::vector::empty
T empty(T... args)
gap_decorator.hpp
Provides seqan3::gap_decorator.
seqan3::alignment_from_cigar
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:81
seqan3::field::alignment
@ alignment
The (pairwise) alignment stored in an object that models seqan3::detail::pairwise_alignment.
Definition record.hpp:78
seqan3::views::slice
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition slice.hpp:137
cigar.hpp
Auxiliary functions for the SAM IO.
std::logic_error
seqan3
The main SeqAn3 namespace.
Definition aligned_sequence_concept.hpp:26
seqan3::get
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:412
std::vector::size
T size(T... args)
slice.hpp
Provides seqan3::views::slice.
std::to_string
T to_string(T... args)
std::tuple
vector
Hide me