io/sam_file/detail/cigar.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2022, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2022, Knut Reinert & MPI für molekulare Genetik
// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
// shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
// -----------------------------------------------------------------------------------------------------
 
#pragma once
 
#include <algorithm>
#include <concepts>
#include <ranges>
#include <seqan3/std/charconv>
#include <sstream>
 
#include <seqan3/alignment/detail/pairwise_alignment_concept.hpp>
#include <seqan3/alphabet/cigar/cigar.hpp>
#include <seqan3/core/debug_stream/detail/to_string.hpp>
#include <seqan3/io/views/detail/take_until_view.hpp>
#include <seqan3/utility/char_operations/predicate.hpp>
#include <seqan3/utility/tuple/concept.hpp>
#include <seqan3/utility/views/single_pass_input.hpp>
#include <seqan3/utility/views/zip.hpp>
 
namespace seqan3::detail
{
struct view_equality_fn
{
    template <std::ranges::forward_range rng1_type, std::ranges::forward_range rng2_type>
    constexpr bool operator()(rng1_type && rng1, rng2_type && rng2) const
    {
        return std::ranges::equal(rng1, rng2);
    }
};
 
inline void update_alignment_lengths(int32_t & ref_length,
                                     int32_t & seq_length,
                                     char const cigar_operation,
                                     uint32_t const cigar_count)
{
    switch (cigar_operation)
    {
    case 'M':
    case '=':
    case 'X':
        ref_length += cigar_count, seq_length += cigar_count;
        break;
    case 'D':
    case 'N':
        ref_length += cigar_count;
        break;
    case 'I':
        seq_length += cigar_count;
        break;
    case 'S':
    case 'H':
    case 'P':
        break; // no op (soft-clipping or padding does not increase either length)
    default:
        throw format_error{"Illegal cigar operation: " + std::string{cigar_operation}};
    }
}
 
template <typename cigar_input_type>
SEQAN3_WORKAROUND_LITERAL std::tuple<std::vector<cigar>, int32_t, int32_t> parse_cigar(cigar_input_type && cigar_input)
{
    std::vector<cigar> operations{};
    std::array<char, 20> buffer{}; // buffer to parse numbers with from_chars. Biggest number should fit in uint64_t
    char cigar_operation{};
    uint32_t cigar_count{};
    int32_t ref_length{}, seq_length{}; // length of aligned part for ref and query
 
    // transform input into a single input view if it isn't already
    auto cigar_view = cigar_input | views::single_pass_input;
 
    // parse the rest of the cigar
    // -------------------------------------------------------------------------------------------------------------
    while (std::ranges::begin(cigar_view) != std::ranges::end(cigar_view)) // until stream is not empty
    {
        auto buff_end = (std::ranges::copy(cigar_view | detail::take_until_or_throw(!is_digit), buffer.data())).out;
        cigar_operation = *std::ranges::begin(cigar_view);
        ++std::ranges::begin(cigar_view);
 
        if (std::from_chars(buffer.begin(), buff_end, cigar_count).ec != std::errc{})
            throw format_error{"Corrupted cigar string encountered"};
 
        update_alignment_lengths(ref_length, seq_length, cigar_operation, cigar_count);
        operations.emplace_back(cigar_count, cigar::operation{}.assign_char(cigar_operation));
    }
 
    return {operations, ref_length, seq_length};
}
 
[[nodiscard]] inline std::string get_cigar_string(std::vector<cigar> const & cigar_vector)
{
    std::string result{};
    std::ranges::for_each(cigar_vector,
                          [&result](auto & cig)
                          {
                              result.append(static_cast<std::string_view>(cig.to_string()));
                          });
    return result;
}
 
template <seqan3::aligned_sequence ref_seq_type, seqan3::aligned_sequence query_seq_type>
[[nodiscard]] inline std::string get_cigar_string(ref_seq_type && ref_seq,
                                                  query_seq_type && query_seq,
                                                  uint32_t const query_start_pos = 0,
                                                  uint32_t const query_end_pos = 0,
                                                  bool const extended_cigar = false)
{
    return get_cigar_string(std::tie(ref_seq, query_seq), query_start_pos, query_end_pos, extended_cigar);
}
 
struct access_restrictor_fn
{
    template <typename chr_t>
    [[noreturn]] chr_t operator()(chr_t) const
    {
        throw std::logic_error{"Access is not allowed because there is no sequence information."};
    }
};
 
} // namespace seqan3::detail