Loading...
Searching...
No Matches
cigar_from_alignment.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/alphabet/cigar/cigar.hpp>
#include <seqan3/alphabet/gap/gap.hpp>
#include <seqan3/utility/views/zip.hpp>
namespace seqan3
{
struct cigar_clipped_bases
{
uint32_t hard_front{};
uint32_t hard_back{};
uint32_t soft_front{};
uint32_t soft_back{};
};
template <typename alignment_type>
inline auto cigar_from_alignment(alignment_type const & alignment,
cigar_clipped_bases const & clipped_bases = {},
bool const extended_cigar = false)
{
static_assert((tuple_like<std::remove_cvref_t<alignment_type>>
&& std::tuple_size_v<std::remove_cvref_t<alignment_type>> == 2),
"The alignment must be a std::pair or std::tuple of size 2.");
static_assert(
(std::equality_comparable_with<gap, std::ranges::range_reference_t<decltype(std::get<0>(alignment))>>
&& std::equality_comparable_with<gap, std::ranges::range_reference_t<decltype(std::get<1>(alignment))>>),
"The alignment must contain two ranges whose value_type is comparable to seqan3::gap.");
/*brief: Compares two seqan3::aligned_sequence values and returns their cigar operation.
* param reference_char The aligned character of the reference to compare.
* param query_char The aligned character of the query to compare.
* param extended_cigar Whether to print the extended cigar alphabet or not. See cigar operation.
* returns A seqan3::cigar::operation representing the alignment operation between the two values.
*
* The resulting CIGAR operation is based on the query character (`query_char`).
*
* ### Example:
*
* The following alignment column shows the reference char ('C') on top and a
* gap for the query char at the bottom.
* ```
* ... C ...
* |
* ... - ...
* ```
* In this case, `to_cigar_op` will return
* 'D' since the query char is "deleted".
*
* The next alignment column shows the reference char ('C') on top and a
* query char ('G') at the bottom.
* ```
* ... C ...
* |
* ... G ...
* ```
* In this case, `to_cigar_op` will return
* 'M', for the basic CIGAR the two bases are aligned, while
* in the extended CIGAR alphabet (`extended_cigar` = `true`) the function
* will return an 'X' since the bases are aligned but are not
* equal.
* \sa seqan3::aligned_sequence
*/
auto to_cigar_op = [extended_cigar](auto const reference_char, auto const query_char)
{
// note that N is not considered because it is equivalent to D but has a special meaning:
// SAM spec: "For mRNA-to-genome alignment, an N operation represents an intron. For other types of alignments,
// the interpretation of N is not defined."
// as we cannot know the meaning, the user has to change D -> N themself
constexpr std::array<char, 6> operators{'M', 'D', 'I', 'P', 'X', '='}; // contains the possible cigar operators.
// no gaps -> 00 -> 0
// gap only in query -> 01 -> 1
// gap only in reference -> 10 -> 2
// both gaps -> 11 -> 3
uint8_t key = (static_cast<uint8_t>(reference_char == gap{}) << 1) | static_cast<uint8_t>(query_char == gap{});
// mismatch -> 100 -> 4
// match -> 101 -> 5
if (extended_cigar && (key == 0)) // in extended format, refine the substitution operator to match/mismatch.
key |= ((1 << 2) | static_cast<uint8_t>(query_char == reference_char)); // maps to [4, 5].
return assign_char_to(operators[key], cigar::operation{});
};
using std::get;
auto & ref_seq = get<0>(alignment);
auto & query_seq = get<1>(alignment);
if (ref_seq.size() != query_seq.size())
throw std::logic_error{"The aligned sequences (including gaps) must have the same length."};
if (std::ranges::empty(ref_seq)) // only check ref_seq because query_seq was checked to have to same size
throw std::logic_error{"The aligned sequences may not be empty."};
std::vector<cigar> result{};
// Add (H)ard-clipping at the start of the query
if (clipped_bases.hard_front)
result.emplace_back(clipped_bases.hard_front, 'H'_cigar_operation);
// Add (S)oft-clipping at the start of the query
if (clipped_bases.soft_front)
result.emplace_back(clipped_bases.soft_front, 'S'_cigar_operation);
// Create cigar string from alignment
// -------------------------------------------------------------------------
// initialize first operation and count value:
cigar::operation operation{to_cigar_op(ref_seq[0], query_seq[0])};
uint32_t count{0};
// go through alignment columns
for (auto && [reference_char, query_char] : views::zip(ref_seq, query_seq))
{
cigar::operation next_op = to_cigar_op(reference_char, query_char);
if (operation == next_op)
{
++count;
}
else
{
result.emplace_back(count, operation);
operation = next_op;
count = 1;
}
}
// append last cigar element
result.emplace_back(count, operation);
// Add (S)oft-clipping at the end of the query
if (clipped_bases.soft_back)
result.emplace_back(clipped_bases.soft_back, 'S'_cigar_operation);
// Add (H)ard-clipping at the end of the query
if (clipped_bases.hard_back)
result.emplace_back(clipped_bases.hard_back, 'H'_cigar_operation);
return result;
}
} // namespace seqan3
cigar.hpp
Provides the seqan3::cigar alphabet.
std::array
seqan3::cigar::operation
exposition_only::cigar_operation operation
The (extended) cigar operation alphabet of M,D,I,H,N,P,S,X,=.
Definition alphabet/cigar/cigar.hpp:93
seqan3::gap
The alphabet of a gap character '-'.
Definition gap.hpp:36
std::vector::emplace_back
T emplace_back(T... args)
gap.hpp
Provides seqan3::gap.
seqan3::assign_char_to
constexpr auto assign_char_to
Assign a character to an alphabet object.
Definition alphabet/concept.hpp:517
seqan3::cigar_from_alignment
auto cigar_from_alignment(alignment_type const &alignment, cigar_clipped_bases const &clipped_bases={}, bool const extended_cigar=false)
Creates a CIGAR string (SAM format) given a seqan3::detail::pairwise_alignment represented by two seq...
Definition cigar_from_alignment.hpp:111
seqan3::field::ref_seq
@ ref_seq
The (reference) "sequence" information, usually a range of nucleotides or amino acids.
Definition record.hpp:80
seqan3::field::alignment
@ alignment
The (pairwise) alignment stored in an object that models seqan3::detail::pairwise_alignment.
Definition record.hpp:78
seqan3::pack_traits::count
constexpr ptrdiff_t count
Count the occurrences of a type in a pack.
Definition type_pack/traits.hpp:161
seqan3::views::zip
seqan::stl::views::zip zip
A view adaptor that takes several views and returns tuple-like values from every i-th element of each...
Definition zip.hpp:24
tuple_like
Whether a type behaves like a tuple.
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
seqan3::cigar_clipped_bases
Helper struct to specialise soft and hard clipping when using seqan3::cigar_from_alignment.
Definition cigar_from_alignment.hpp:27
seqan3::cigar_clipped_bases::hard_front
uint32_t hard_front
The number of hard clipped bases at the front of the CIGAR string.
Definition cigar_from_alignment.hpp:28
seqan3::cigar_clipped_bases::soft_back
uint32_t soft_back
The number of soft clipped bases at the back of the CIGAR string.
Definition cigar_from_alignment.hpp:31
seqan3::cigar_clipped_bases::hard_back
uint32_t hard_back
The number of hard clipped bases at the back of the CIGAR string.
Definition cigar_from_alignment.hpp:29
seqan3::cigar_clipped_bases::soft_front
uint32_t soft_front
The number of soft clipped bases at the front of the CIGAR string.
Definition cigar_from_alignment.hpp:30
std::tuple_size_v
T tuple_size_v
vector
zip.hpp
Provides seqan3::views::zip.
Hide me