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}};
    }
}
 
SEQAN3_WORKAROUND_LITERAL std::vector<cigar> parse_cigar(std::string_view const cigar_str)
{
    std::vector<seqan3::cigar> cigar_vector{};
 
    if (cigar_str == "*")
        return cigar_vector;
 
    uint32_t cigar_count{};
    char const * ptr = cigar_str.data();
    char const * const end = ptr + cigar_str.size();
 
    while (ptr < end)
    {
        auto const res = std::from_chars(ptr, end, cigar_count); // reads number up to next character
 
        if (res.ec != std::errc{})
            throw format_error{"Corrupted cigar string."};
 
        ptr = res.ptr + 1; // skip cigar operation character
 
        cigar_vector.emplace_back(cigar_count, seqan3::assign_char_strictly_to(*res.ptr, seqan3::cigar::operation{}));
    }
 
    return cigar_vector;
}
 
[[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