detail.hpp

Go to the documentation of this file.

// -----------------------------------------------------------------------------------------------------
// Copyright (c) 2006-2019, Knut Reinert & Freie Universität Berlin
// Copyright (c) 2016-2019, 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 <sstream>

#include <seqan3/alignment/aligned_sequence/aligned_sequence_concept.hpp>
#include <seqan3/core/concept/tuple.hpp>
#include <seqan3/range/view/single_pass_input.hpp>
#include <seqan3/range/view/take_until.hpp>
#include <seqan3/std/algorithm>
#include <seqan3/std/charconv>
#include <seqan3/std/concepts>
#include <seqan3/std/ranges>

#include <range/v3/view/zip.hpp>

namespace seqan3::detail
{
struct view_equality_fn
{
    template <std::ranges::ForwardRange rng1_type, std::ranges::ForwardRange rng2_type>
    constexpr bool operator()(rng1_type && rng1, rng2_type && rng2) const
    {
        return std::ranges::equal(rng1, rng2);
    }
};

template <typename reference_char_type, typename query_char_type>
    requires std::detail::WeaklyEqualityComparableWith<reference_char_type, gap> &&
             std::detail::WeaklyEqualityComparableWith<query_char_type, gap>
char compare_aligned_values(reference_char_type const reference_char,
                            query_char_type const query_char,
                            bool const extended_cigar)
{
    return (reference_char == gap{})
                ? (query_char == gap{})
                    ? 'P'
                    : 'I'
                : (query_char == gap{})
                    ? 'D'
                    : (extended_cigar)
                        ? (query_char == reference_char)
                            ? '='
                            : 'X'
                        : 'M';
}

template <std::ranges::ForwardRange ref_seq_type, std::ranges::ForwardRange query_seq_type>
    requires std::detail::WeaklyEqualityComparableWith<gap, reference_t<ref_seq_type>> &&
             std::detail::WeaklyEqualityComparableWith<gap, reference_t<query_seq_type>>
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)
{
    if (ref_seq.size() != query_seq.size())
        throw std::logic_error{"The aligned sequences must have the same length."};

    std::ostringstream result;

    if (!ref_seq.size())
        return std::string(); // return empty string if sequences are empty

    // Add (S)oft-clipping at the start of the read
    if (query_start_pos)
        result << query_start_pos << 'S';

    // Create cigar string from alignment
    // -------------------------------------------------------------------------
    // initialize first operation:
    char tmp_char{compare_aligned_values(ref_seq[0], query_seq[0], extended_cigar)};
    size_t tmp_length{0};

    // go through alignment columns
    for (auto column : std::view::zip(ref_seq, query_seq))
    {
        char next_op = compare_aligned_values(std::get<0>(column), std::get<1>(column), extended_cigar);

        if (tmp_char == next_op)
        {
            ++tmp_length;
        }
        else
        {
            result << tmp_length << tmp_char;
            tmp_char = next_op;
            tmp_length = 1;
        }
    }
    // append last cigar element
    result << tmp_length << tmp_char;

    // Add (S)oft-clipping at the end of the read
    if (query_end_pos)
        result << query_end_pos << 'S';

    return result.str();
}

template <TupleLike alignment_type>
    requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2
std::string get_cigar_string(alignment_type && alignment,
                             uint32_t const query_start_pos = 0,
                             uint32_t const query_end_pos = 0,
                             bool const extended_cigar = false)
{
    return get_cigar_string(get<0>(alignment), get<1>(alignment),
                            query_start_pos, query_end_pos, extended_cigar);
}
template <TupleLike alignment_type>
    requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2
std::vector<std::pair<char, size_t>> get_cigar_vector(alignment_type && alignment,
                                                      uint32_t const query_start_pos = 0,
                                                      uint32_t const query_end_pos = 0,
                                                      bool const extended_cigar = false)
{
    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 must have the same length."};

    std::vector<std::pair<char, size_t>> result{};

    if (!ref_seq.size())
        return result; // return empty string if sequences are empty

    // Add (S)oft-clipping at the start of the read
    if (query_start_pos)
        result.emplace_back('S', query_start_pos);

    // Create cigar string from alignment
    // -------------------------------------------------------------------------
    // initialize first operation:
    char tmp_char{compare_aligned_values(ref_seq[0], query_seq[0], extended_cigar)};
    size_t tmp_length{0};

    // go through alignment columns
    for (auto column : std::view::zip(ref_seq, query_seq))
    {
        char next_op = compare_aligned_values(std::get<0>(column), std::get<1>(column), extended_cigar);

        if (tmp_char == next_op)
        {
            ++tmp_length;
        }
        else
        {
            result.emplace_back(tmp_char, tmp_length);
            tmp_char = next_op;
            tmp_length = 1;
        }
    }

    // append last cigar element
    result.emplace_back(tmp_char, tmp_length);

    // Add (S)oft-clipping at the end of the read
    if (query_end_pos)
        result.emplace_back('S', query_end_pos);

    return result;
}

template <std::ranges::InputRange cigar_input_type>
std::tuple<std::vector<std::pair<char, size_t>>, size_t, size_t, size_t, size_t>
parse_cigar(cigar_input_type && cigar_input)
{
    std::vector<std::pair<char, size_t>> operations{};
    size_t sc_begin_count{}; // number of soft clipped bases at the beginning
    size_t sc_end_count{};   // number of soft clipped bases at the end
    std::array<char, 20> buffer{}; // buffer to parse numbers with from_chars. Biggest number should fit in uint64_t
    char cigar_op{'\0'};
    size_t cigar_count{}, ref_length{}, seq_length{}; // length of aligned part for ref and query

    auto update_lengths_fn = [&ref_length, &seq_length, &cigar_op, &cigar_count] ()
        {
            if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
            {
                ref_length += cigar_count;
                seq_length += cigar_count;
            }
            else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op))
            {
                ref_length += cigar_count;
            }
            else if (is_char<'I'>(cigar_op))
            {
                seq_length += cigar_count;
            }
            else // illegal character
            {
                if (is_char<'P'>(cigar_op))
                    throw format_error{"We do currently not support cigar operation 'P'."};
                else
                    throw format_error{std::string{"Illegal cigar operation: "} + std::string{cigar_op}};
            }
        };

    // transform input into a signle input view if it isn't already
    auto cigar_view = cigar_input | view::single_pass_input;

    // check hard/soft clipping at the beginning manually
    // -----------------------------------------------------------------------------------------------------------------
    auto [ignore, buffer_end] = std::ranges::copy(cigar_view | view::take_until_or_throw(!is_digit), buffer.data());
    (void) ignore;

    cigar_op = *std::ranges::begin(cigar_view);
    std::ranges::next(std::ranges::begin(cigar_view));

    if (is_char<'H'>(cigar_op)) // hard clipping is ignored. parse the next operation
    {
        auto [ignore2, buffer_end2] = std::ranges::copy(cigar_view
                                                        | view::take_until_or_throw(!is_digit), buffer.data());
        buffer_end = buffer_end2;
        (void) ignore2;

        cigar_op = *std::ranges::begin(cigar_view);
        std::ranges::next(std::ranges::begin(cigar_view));
    }

    if (std::from_chars(buffer.begin(), buffer_end, cigar_count).ec != std::errc{})
        throw format_error{"Corrupted cigar string encountered"};

    if (is_char<'S'>(cigar_op)) // check for soft clipping at the beginning
    {
        sc_begin_count = cigar_count;
    }
    else
    {
        update_lengths_fn();
        operations.push_back({cigar_op, cigar_count});
    }

    // parse the rest of the cigar
    // -----------------------------------------------------------------------------------------------------------------
    while (std::ranges::begin(cigar_view) != std::ranges::end(cigar_view)) // until stream is not empty
    {
        buffer_end = (std::ranges::copy(cigar_view | view::take_until_or_throw(!is_digit), buffer.data())).out;
        cigar_op = *std::ranges::begin(cigar_view);
        std::ranges::next(std::ranges::begin(cigar_view));

        if (std::from_chars(buffer.begin(), buffer_end, cigar_count).ec != std::errc{})
            throw format_error{"Corrupted cigar string encountered"};

        if (is_char<'S'>(cigar_op)) // we are at the end, hard clipping afterwards can be ignored
        {
            sc_end_count = cigar_count;
            return {operations, ref_length, seq_length, sc_begin_count, sc_end_count};
        }
        update_lengths_fn();
        operations.push_back({cigar_op, cigar_count});
    }
    return {operations, ref_length, seq_length, sc_begin_count, sc_end_count};
}

template <std::ranges::InputRange cigar_input_type>
auto parse_binary_cigar(cigar_input_type && cigar_input, uint16_t n_cigar_op)
{
    std::vector<std::pair<char, size_t>> operations{};
    size_t sc_begin_count{}; // number of soft clipped bases at the beginning
    size_t sc_end_count{};   // number of soft clipped bases at the end
    char cigar_op{'\0'};
    size_t cigar_count{},ref_length{}, seq_length{};
    uint32_t o_and_c{};
    static constexpr char const * CIGAR_MAPPING = "MIDNSHP=X*******";
    static constexpr uint32_t CIGAR_MASK = 0x0f; // 0000000000001111

    auto update_lengths_fn = [&ref_length, &seq_length, &cigar_op, &cigar_count] ()
        {
            if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
            {
                ref_length += cigar_count;
                seq_length += cigar_count;
            }
            else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op))
            {
                ref_length += cigar_count;
            }
            else if (is_char<'I'>(cigar_op))
            {
                seq_length += cigar_count;
            }
            else // illegal character
            {
                if (is_char<'P'>(cigar_op))
                    throw format_error{"We do currently not support cigar operation 'P'."};
                else
                    throw format_error{std::string{"Illegal cigar operation: "} + std::string{cigar_op}};
            }
        };

    if (n_cigar_op == 0) // [[unlikely]]
        return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};

    std::ranges::copy_n(std::ranges::begin(cigar_input), sizeof(o_and_c), reinterpret_cast<char*>(&o_and_c));
    assert((o_and_c & CIGAR_MASK) <= 8u);
    cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
    cigar_count = o_and_c >> 4;
    --n_cigar_op;

    if (is_char<'H'>(cigar_op)) // hard clipping is ignored. parse the next operation
    {
        if (n_cigar_op == 0)
            return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};

        std::ranges::copy_n(std::ranges::begin(cigar_input), sizeof(o_and_c), reinterpret_cast<char*>(&o_and_c));
        assert((o_and_c & CIGAR_MASK) <= 8u);
        cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
        cigar_count = o_and_c >> 4;
        --n_cigar_op;
    }

    if (is_char<'S'>(cigar_op)) // check for soft clipping at the beginning
    {
        sc_begin_count = cigar_count;
    }
    else
    {
        update_lengths_fn();
        operations.push_back({cigar_op, cigar_count});
    }

    // parse the rest of the cigar
    // -----------------------------------------------------------------------------------------------------------------
    while (n_cigar_op > 0) // until stream is not empty
    {
        std::ranges::copy_n(std::ranges::begin(cigar_input), sizeof(o_and_c), reinterpret_cast<char*>(&o_and_c));
        assert((o_and_c & CIGAR_MASK) <= 8u);
        cigar_op = CIGAR_MAPPING[o_and_c & CIGAR_MASK];
        cigar_count = o_and_c >> 4;

        if (is_char<'S'>(cigar_op)) // we are at the end, hard clipping afterwards can be ignored
        {
            sc_end_count = cigar_count;
            return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
        }

        update_lengths_fn();
        operations.push_back({cigar_op, cigar_count});
        --n_cigar_op;
    }
    return std::tuple{operations, ref_length, seq_length, sc_begin_count, sc_end_count};
}

template <TupleLike alignment_type>
    requires std::tuple_size_v<remove_cvref_t<alignment_type>> == 2 &&
             detail::all_satisfy_aligned_seq<detail::tuple_type_list_t<alignment_type>>
void alignment_from_cigar(alignment_type & alignment, std::vector<std::pair<char, size_t>> const & cigar)
{
    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_op, cigar_count] : cigar)
    {
        if (is_char<'M'>(cigar_op) || is_char<'='>(cigar_op) || is_char<'X'>(cigar_op))
        {
            std::ranges::advance(current_ref_pos , cigar_count);
            std::ranges::advance(current_read_pos, cigar_count);
        }
        else if (is_char<'D'>(cigar_op) || is_char<'N'>(cigar_op)) // insert gaps into read
        {
            assert(std::distance(current_read_pos, std::ranges::end(get<1>(alignment))) >= 0);
            current_read_pos = insert_gap(get<1>(alignment), current_read_pos, cigar_count);
            ++current_read_pos;
            std::ranges::advance(current_ref_pos , cigar_count);
        }
        else if (is_char<'I'>(cigar_op)) // Insert gaps into ref
        {
            assert(std::ranges::distance(current_ref_pos, std::ranges::end(get<0>(alignment))) >= 0);
            current_ref_pos = insert_gap(get<0>(alignment), current_ref_pos, cigar_count);
            ++current_ref_pos;
            std::ranges::advance(current_read_pos, cigar_count);
        }
        else // illegal character
        {
            if (is_char<'P'>(cigar_op))
                throw format_error{"We do currently not support cigar operation 'P'."};
            else
                throw format_error{std::string{"Illegal cigar operation: "} + std::string{cigar_op}};
        }
    }
}

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