aligned_sequence_concept.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 <iomanip>
#include <tuple>

#include <range/v3/algorithm/for_each.hpp>

#include <seqan3/alignment/exception.hpp>
#include <seqan3/alphabet/gap/all.hpp>
#include <seqan3/core/concept/tuple.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/core/type_traits/all.hpp>
#include <seqan3/range/container/concept.hpp>
#include <seqan3/range/view/slice.hpp>
#include <seqan3/range/view/to_char.hpp>
#include <seqan3/std/ranges>

// ---------------------------------------------------------------------------------------------------------------------
// unaligned_seq transformation trait
// ---------------------------------------------------------------------------------------------------------------------

namespace seqan3::detail
{

template <template <typename...> typename container_type, typename seq_alph_t, typename ...rest_t>
    requires Container<container_type<gapped<seq_alph_t>, rest_t...>>
constexpr auto remove_gap_from_value_type(container_type<gapped<seq_alph_t>, rest_t...>)
    -> container_type<seq_alph_t, rest_t...>;

template <template <typename...> typename container_type,
          template <typename...> typename allocator_type,
          typename seq_alph_t, typename ...rest_t>
    requires Container<container_type<gapped<seq_alph_t>, allocator_type<gapped<seq_alph_t>>, rest_t...>>
constexpr auto remove_gap_from_value_type(container_type<gapped<seq_alph_t>, allocator_type<gapped<seq_alph_t>>, rest_t...>)
    -> container_type<seq_alph_t, allocator_type<seq_alph_t>, rest_t...>;

template <typename t>
struct unaligned_seq
{};

template <typename t>
    requires !requires { typename std::remove_reference_t<t>::unaligned_seq_type; } &&
              requires { remove_gap_from_value_type(std::declval<t>()); }
struct unaligned_seq<t>
{
    using type = decltype(remove_gap_from_value_type(std::declval<t>()));
};

// customisation point for our gap decorators.
template <typename t>
    requires requires { typename std::remove_reference_t<t>::unaligned_seq_type; }
struct unaligned_seq<t>
{
    using type = typename std::remove_reference_t<t>::unaligned_seq_type; 
};

template <typename t>
using unaligned_seq_t = typename unaligned_seq<t>::type;

} // namespace seqan3::detail

// ---------------------------------------------------------------------------------------------------------------------
// AlignedSequence
// ---------------------------------------------------------------------------------------------------------------------

namespace seqan3
{

template <typename t>
SEQAN3_CONCEPT AlignedSequence =
    std::ranges::ForwardRange<t> &&
    Alphabet<reference_t<t>> &&
    WeaklyAssignable<reference_t<t>, gap const &> &&
    requires { typename detail::unaligned_seq_t<t>; } &&
    requires (t v, detail::unaligned_seq_t<t> unaligned)
    {
        { insert_gap(v, v.begin()) } -> typename t::iterator; // global functions for generic usability
        { insert_gap(v, v.begin(), 2) } -> typename t::iterator;
        { erase_gap(v, v.begin()) } -> typename t::iterator;
        { erase_gap(v, v.begin(), v.end()) } -> typename t::iterator;
        { assign_unaligned(v, unaligned) } -> void;
    };

// ---------------------------------------------------------------------------------------------------------------------
// Aligned sequence interface for containers over the seqan3::gapped alphabet
// ---------------------------------------------------------------------------------------------------------------------

template <SequenceContainer aligned_seq_t>
    requires detail::is_gapped_alphabet<value_type_t<aligned_seq_t>>
inline typename aligned_seq_t::iterator insert_gap(aligned_seq_t & aligned_seq,
                                                   typename aligned_seq_t::const_iterator pos_it)
{
    return aligned_seq.insert(pos_it, value_type_t<aligned_seq_t>{gap{}});
}

template <SequenceContainer aligned_seq_t>
    requires detail::is_gapped_alphabet<value_type_t<aligned_seq_t>>
inline typename aligned_seq_t::iterator insert_gap(aligned_seq_t & aligned_seq,
                                                   typename aligned_seq_t::const_iterator pos_it,
                                                   typename aligned_seq_t::size_type size)
{
    return aligned_seq.insert(pos_it, size, value_type_t<aligned_seq_t>{gap{}});
}

template <SequenceContainer aligned_seq_t>
    requires detail::is_gapped_alphabet<value_type_t<aligned_seq_t>>
inline typename aligned_seq_t::iterator erase_gap(aligned_seq_t & aligned_seq,
                                                  typename aligned_seq_t::const_iterator pos_it)
{
    if (*pos_it != gap{}) // [[unlikely]]
        throw gap_erase_failure("The position to be erased does not contain a gap.");

    return aligned_seq.erase(pos_it);
}

template <SequenceContainer aligned_seq_t>
    requires detail::is_gapped_alphabet<value_type_t<aligned_seq_t>>
inline typename aligned_seq_t::iterator erase_gap(aligned_seq_t & aligned_seq,
                                                  typename aligned_seq_t::const_iterator first,
                                                  typename aligned_seq_t::const_iterator last)
{
    for (auto it = first; it != last; ++it)
        if (*it != gap{}) // [[unlikely]]
            throw gap_erase_failure("The range to be erased contains at least one non-gap character.");

    return aligned_seq.erase(first, last);
}

template <SequenceContainer aligned_seq_t, std::ranges::ForwardRange unaligned_sequence_type>
    requires detail::is_gapped_alphabet<value_type_t<aligned_seq_t>> &&
             WeaklyAssignable<reference_t<aligned_seq_t>, reference_t<unaligned_sequence_type>>
inline void assign_unaligned(aligned_seq_t & aligned_seq, unaligned_sequence_type && unaligned_seq)
{
    using std::swap;
    aligned_seq_t tmp;
    tmp.resize(std::ranges::distance(unaligned_seq));
    std::copy(std::ranges::begin(unaligned_seq), std::ranges::end(unaligned_seq), std::ranges::begin(tmp));
    swap(aligned_seq, tmp);
}

template <typename range_type>
    requires requires (range_type v)
        {
            v.insert_gap(typename range_type::iterator{});
            v.insert_gap(typename range_type::iterator{}, typename range_type::size_type{});
        }
typename range_type::iterator insert_gap(range_type & rng,
                                         typename range_type::iterator const it,
                                         typename range_type::size_type const size = 1)
{
    return rng.insert_gap(it, size);
}

template <typename range_type>
    requires requires (range_type v) { v.erase_gap(typename range_type::iterator{}); }
typename range_type::iterator erase_gap(range_type & rng,
                                        typename range_type::iterator const it)
{
    return rng.erase_gap(it);
}

template <typename range_type>
    requires requires (range_type v) { v.erase_gap(typename range_type::iterator{}, typename range_type::iterator{}); }
typename range_type::iterator erase_gap(range_type & rng,
                                        typename range_type::iterator const first,
                                        typename range_type::iterator const last)
{
    return rng.erase_gap(first, last);
}

namespace detail
{

template<TupleLike alignment_t, size_t ...idx>
void stream_alignment(debug_stream_type & stream, alignment_t const & align, std::index_sequence<idx...> const & )
{
    using std::get;
    size_t const alignment_size = get<0>(align).size();

    // split alignment into blocks of length 50 and loop over parts
    for (size_t begin_pos = 0; begin_pos < alignment_size; begin_pos += 50)
    {
        size_t const end_pos = std::min(begin_pos + 50, alignment_size);

        // write header line
        if (begin_pos != 0)
            stream << '\n';

        stream << std::setw(7) << begin_pos << ' ';
        for (size_t pos = begin_pos + 1; pos <= end_pos; ++pos)
        {
            if (pos % 10 == 0)
                stream << ':';
            else if (pos % 5 == 0)
                stream << '.';
            else
                stream << ' ';
        }

        // write first sequence
        stream << '\n' << std::setw(8) << "";
        ranges::for_each(get<0>(align) | view::slice(begin_pos, end_pos) | view::to_char,
                         [&stream] (char ch) { stream << ch; });

        auto stream_f = [&] (auto const & previous_seq, auto const & aligned_seq)
        {
            // write alignment bars
            stream << '\n' << std::setw(8) << "";
            ranges::for_each(ranges::zip_view(previous_seq, aligned_seq) | view::slice(begin_pos, end_pos),
                             [&stream] (auto && ch) { stream << (get<0>(ch) == get<1>(ch) ? '|' : ' '); });

            // write next sequence
            stream << '\n' << std::setw(8) << "";
            ranges::for_each(aligned_seq | view::slice(begin_pos, end_pos) | view::to_char,
                             [&stream] (char ch) { stream << ch; });
        };
        (stream_f(get<idx>(align), get<idx + 1>(align)), ...);
        stream << '\n';
    }
}

template <typename ...elems>
inline bool constexpr all_satisfy_aligned_seq = false;

template <typename ...elems>
inline bool constexpr all_satisfy_aligned_seq<type_list<elems...>> = (AlignedSequence<elems> && ...);

} // namespace detail

template <TupleLike tuple_t>
    requires detail::all_satisfy_aligned_seq<detail::tuple_type_list_t<tuple_t>>
inline debug_stream_type & operator<<(debug_stream_type & stream, tuple_t const & alignment)
{
    static_assert(std::tuple_size_v<tuple_t> >= 2, "An alignment requires at least two sequences.");
    detail::stream_alignment(stream, alignment, std::make_index_sequence<std::tuple_size_v<tuple_t> - 1> {});
    return stream;
}

} // namespace seqan