format_sam.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 <iterator>
#include <string>
#include <vector>

#include <seqan3/core/char_operations/predicate.hpp>
#include <seqan3/core/concept/core_language.hpp>
#include <seqan3/core/concept/tuple.hpp>
#include <seqan3/core/detail/reflection.hpp>
#include <seqan3/core/detail/to_string.hpp>
#include <seqan3/core/type_traits/range.hpp>
#include <seqan3/core/type_traits/template_inspection.hpp>
#include <seqan3/io/alignment_file/detail.hpp>
#include <seqan3/io/alignment_file/header.hpp>
#include <seqan3/io/alignment_file/input_format_concept.hpp>
#include <seqan3/io/alignment_file/input_options.hpp>
#include <seqan3/io/alignment_file/output_format_concept.hpp>
#include <seqan3/io/alignment_file/output_options.hpp>
#include <seqan3/io/alignment_file/sam_tag_dictionary.hpp>
#include <seqan3/io/detail/ignore_output_iterator.hpp>
#include <seqan3/io/detail/misc.hpp>
#include <seqan3/io/stream/iterator.hpp>
#include <seqan3/range/decorator/gap_decorator.hpp>
#include <seqan3/range/detail/misc.hpp>
#include <seqan3/range/view/char_to.hpp>
#include <seqan3/range/view/istreambuf.hpp>
#include <seqan3/range/view/repeat_n.hpp>
#include <seqan3/range/view/slice.hpp>
#include <seqan3/range/view/take_until.hpp>
#include <seqan3/range/view/to_char.hpp>
#include <seqan3/std/algorithm>
#include <seqan3/std/charconv>
#include <seqan3/std/concepts>
#include <seqan3/std/ranges>
#include <seqan3/version.hpp>

namespace seqan3
{

struct format_sam
{
    static inline std::vector<std::string> file_extensions
    {
        { "sam" },
    };

    static constexpr char format_version[4] = "1.6";
};

} // namespace seqan3

namespace seqan3::detail
{

template <>
class alignment_file_input_format<format_sam>
{
public:
    using format_tag = format_sam;

    alignment_file_input_format()                                                noexcept = default; 
    alignment_file_input_format(alignment_file_input_format const &)                      = delete;
    alignment_file_input_format & operator=(alignment_file_input_format const &)          = delete;
    alignment_file_input_format(alignment_file_input_format &&)                  noexcept = default; 
    alignment_file_input_format & operator=(alignment_file_input_format &&)      noexcept = default; 
    ~alignment_file_input_format()                                               noexcept = default; 

    template <typename stream_type,     // constraints checked by file
              typename seq_legal_alph_type,
              typename ref_seqs_type,
              typename ref_ids_type,
              typename seq_type,
              typename id_type,
              typename offset_type,
              typename ref_seq_type,
              typename ref_id_type,
              typename ref_offset_type,
              typename align_type,
              typename flag_type,
              typename mapq_type,
              typename qual_type,
              typename mate_type,
              typename tag_dict_type,
              typename e_value_type,
              typename bit_score_type>
    void read(stream_type                                             & stream,
              alignment_file_input_options<seq_legal_alph_type> const & SEQAN3_DOXYGEN_ONLY(options),
              ref_seqs_type                                           & ref_seqs,
              alignment_file_header<ref_ids_type>                     & header,
              seq_type                                                & seq,
              qual_type                                               & qual,
              id_type                                                 & id,
              offset_type                                             & offset,
              ref_seq_type                                            & SEQAN3_DOXYGEN_ONLY(ref_seq),
              ref_id_type                                             & ref_id,
              ref_offset_type                                         & ref_offset,
              align_type                                              & align,
              flag_type                                               & flag,
              mapq_type                                               & mapq,
              mate_type                                               & mate,
              tag_dict_type                                           & tag_dict,
              e_value_type                                            & SEQAN3_DOXYGEN_ONLY(e_value),
              bit_score_type                                          & SEQAN3_DOXYGEN_ONLY(bit_score))
    {
        static_assert(detail::decays_to_ignore_v<ref_offset_type> ||
                      detail::is_type_specialisation_of_v<ref_offset_type, std::optional>,
                      "The ref_offset must be a specialisation of std::optional.");

        auto stream_view = view::istreambuf(stream);
        auto field_view = stream_view | view::take_until_or_throw_and_consume(is_char<'\t'>);

        // these variables need to be stored to compute the ALIGNMENT
        int32_t ref_offset_tmp{};
        value_type_t<decltype(header.ref_ids())> ref_id_tmp{};
        [[maybe_unused]] int32_t offset_tmp{};
        [[maybe_unused]] int32_t soft_clipping_end{};
        [[maybe_unused]] std::vector<std::pair<char, size_t>> cigar{};
        [[maybe_unused]] int32_t ref_length{0}, seq_length{0}; // length of aligned part for ref and query

        // Header
        // -------------------------------------------------------------------------------------------------------------
        if (is_char<'@'>(*std::ranges::begin(stream_view))) // we always read the header if present
        {
            read_header(stream_view, header, ref_seqs);

            if (std::ranges::begin(stream_view) == std::ranges::end(stream_view)) // file has no records
                return;
        }

        // Fields 1-5: ID FLAG REF_ID REF_OFFSET MAPQ
        // -------------------------------------------------------------------------------------------------------------
        read_field(field_view, id);

        read_field(field_view, flag);

        read_field(field_view, ref_id_tmp);
        check_and_assign_ref_id(ref_id, ref_id_tmp, header, ref_seqs);

        read_field(field_view, ref_offset_tmp);
        --ref_offset_tmp; // SAM format is 1-based but SeqAn operates 0-based

        if (ref_offset_tmp == -1)
            ref_offset = std::nullopt; // indicates an unmapped read -> ref_offset is not set
        else if (ref_offset_tmp > -1)
            ref_offset = ref_offset_tmp;
        else if (ref_offset_tmp < -1)
            throw format_error{"No negative values are allowed for field::REF_OFFSET."};

        read_field(field_view, mapq);

        // Field 6: CIGAR
        // -------------------------------------------------------------------------------------------------------------
        if constexpr (!detail::decays_to_ignore_v<align_type>)
        {
            if (!is_char<'*'>(*std::ranges::begin(stream_view))) // no cigar information given
            {
                std::tie(cigar, ref_length, seq_length, offset_tmp, soft_clipping_end) = parse_cigar(field_view);
            }
            else
            {
                std::ranges::next(std::ranges::begin(field_view)); // skip '*'
            }
        }
        else
        {
            detail::consume(field_view);
        }

        offset = offset_tmp;

        // Field 7-9: (RNEXT PNEXT TLEN) = MATE
        // -------------------------------------------------------------------------------------------------------------
        if constexpr (!detail::decays_to_ignore_v<mate_type>)
        {
            value_type_t<decltype(header.ref_ids())> tmp_mate_ref_id{};
            read_field(field_view, tmp_mate_ref_id); // RNEXT

            if (tmp_mate_ref_id == "=") // indicates "same as ref id"
            {
                if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
                    get<0>(mate) = ref_id;
                else
                    check_and_assign_ref_id(get<0>(mate), ref_id_tmp, header, ref_seqs);
            }
            else
            {
                check_and_assign_ref_id(get<0>(mate), tmp_mate_ref_id, header, ref_seqs);
            }

            int32_t tmp_pnext{};
            read_field(field_view, tmp_pnext); // PNEXT

            if (tmp_pnext > 0)
                get<1>(mate) = --tmp_pnext; // SAM format is 1-based but SeqAn operates 0-based.
            else if (tmp_pnext < 0)
                throw format_error{"No negative values are allowed at the mate mapping position."};
            // tmp_pnext == 0 indicates an unmapped mate -> do not fill std::optional get<1>(mate)

            read_field(field_view, get<2>(mate)); // TLEN
        }
        else
        {
            for (size_t i = 0; i < 3u; ++i)
            {
                detail::consume(field_view);
            }
        }

        // Field 10: Sequence
        // -------------------------------------------------------------------------------------------------------------
        if (!is_char<'*'>(*std::ranges::begin(stream_view))) // sequence information is given
        {
            auto constexpr is_legal_alph = is_in_alphabet<seq_legal_alph_type>;
            auto seq_stream = field_view | std::view::transform([is_legal_alph] (char const c) // enforce legal alphabet
                                           {
                                               if (!is_legal_alph(c))
                                                   throw format_error{std::string{"Encountered an unexpected letter: "} +
                                                                     is_legal_alph.msg.str() +
                                                                     " evaluated to false on " +
                                                                     detail::make_printable(c)};
                                               return c;
                                           });

            if constexpr (detail::decays_to_ignore_v<seq_type>)
            {
                if constexpr (!detail::decays_to_ignore_v<align_type>)
                {
                    static_assert(SequenceContainer<std::remove_reference_t<decltype(get<1>(align))>>,
                                  "If you want to read ALIGNMENT but not SEQ, the alignment"
                                  " object must store a sequence container at the second (query) position.");

                    if (!cigar.empty()) // only parse alignment if cigar information was given
                    {

                        auto tmp_iter = std::ranges::begin(seq_stream);
                        std::ranges::advance(tmp_iter, offset_tmp);

                        for (; seq_length > 0; --seq_length) // seq_length is not needed anymore
                        {
                            get<1>(align).push_back(value_type_t<decltype(get<1>(align))>{}.assign_char(*tmp_iter));
                            ++tmp_iter;
                        }

                        std::ranges::advance(tmp_iter, soft_clipping_end);
                    }
                    else
                    {
                        get<1>(align) = std::remove_reference_t<decltype(get<1>(align))>{}; // empty container
                    }
                }
                else
                {
                    detail::consume(seq_stream);
                }
            }
            else
            {
                read_field(seq_stream, seq);

                if constexpr (!detail::decays_to_ignore_v<align_type>)
                {
                    if (!cigar.empty()) // if no alignment info is given, the field::ALIGNMENT should remain empty
                    {
                        assign_unaligned(get<1>(align),
                                         seq | view::slice(static_cast<decltype(std::ranges::size(seq))>(offset_tmp),
                                                           std::ranges::size(seq) - soft_clipping_end));
                    }
                }
            }
        }
        else
        {
            std::ranges::next(std::ranges::begin(field_view)); // skip '*'
        }

        // Field 11:  Quality
        // -------------------------------------------------------------------------------------------------------------
        auto const tab_or_end = is_char<'\t'> || is_char<'\r'> || is_char<'\n'>;
        read_field(stream_view | view::take_until_or_throw(tab_or_end), qual);

        if constexpr (!detail::decays_to_ignore_v<seq_type> && !detail::decays_to_ignore_v<qual_type>)
        {
            if (std::ranges::distance(seq) != 0 && std::ranges::distance(qual) != 0 &&
                std::ranges::distance(seq) != std::ranges::distance(qual))
            {
                throw format_error{to_string("Sequence length (", std::ranges::distance(seq), ") and quality length (",
                                             std::ranges::distance(qual), ") must be the same.")};
            }
        }

        // All remaining optional fields if any: SAM tags dictionary
        // -------------------------------------------------------------------------------------------------------------
        while (is_char<'\t'>(*std::ranges::begin(stream_view))) // read all tags if present
        {
            std::ranges::next(std::ranges::begin(stream_view)); // skip tab
            read_field(stream_view | view::take_until_or_throw(tab_or_end), tag_dict);
        }

        detail::consume(stream_view | view::take_until(!(is_char<'\r'> || is_char<'\n'>))); // consume new line

        // DONE READING - wrap up
        // -------------------------------------------------------------------------------------------------------------
        // Alignment object construction
        // Note that the query sequence in get<1>(align) has already been filled while reading Field 10.
        if constexpr (!detail::decays_to_ignore_v<align_type>)
        {
            int32_t ref_idx{(ref_id_tmp.empty()/*unmapped read?*/) ? -1 : 0};

            if constexpr (!detail::decays_to_ignore_v<ref_seqs_type>)
            {
                if (!ref_id_tmp.empty())
                {
                    assert(header.ref_dict.count(ref_id_tmp) != 0); // taken care of in check_and_assign_ref_id()
                    ref_idx = header.ref_dict[ref_id_tmp];          // get index for reference sequence
                }
            }

            construct_alignment(align, cigar, ref_idx, ref_seqs, ref_offset_tmp, ref_length);
        }
    }

protected:
    std::array<char, 316> buffer{}; // Doubles can be up to 316 characters

    bool ref_info_present_in_header{false};

    template <typename ref_id_type,
              typename ref_id_tmp_type,
              typename header_type,
              typename ref_seqs_type>
    void check_and_assign_ref_id(ref_id_type      & ref_id,
                                 ref_id_tmp_type  & ref_id_tmp,
                                 header_type      & header,
                                 ref_seqs_type    & /*tag*/)
    {
        if (!std::ranges::empty(ref_id_tmp)) // otherwise the std::optional will not be filled
        {
            auto search = header.ref_dict.find(ref_id_tmp);

            if (search == header.ref_dict.end())
            {
                if constexpr(detail::decays_to_ignore_v<ref_seqs_type>) // no reference information given
                {
                    if (ref_info_present_in_header)
                    {
                        throw format_error{"Unknown reference id found in record which is not present in the header."};
                    }
                    else
                    {
                        header.ref_ids().push_back(ref_id_tmp);
                        auto pos = std::ranges::size(header.ref_ids()) - 1;
                        header.ref_dict[header.ref_ids()[pos]] = pos;
                        ref_id = pos;
                    }
                }
                else
                {
                    throw format_error{"Unknown reference id found in record which is not present in the given ids."};
                }
            }
            else
            {
                ref_id = search->second;
            }
        }
    }

    template <typename align_type, typename ref_seqs_type>
    void construct_alignment(align_type                           & align,
                             std::vector<std::pair<char, size_t>> & cigar,
                             [[maybe_unused]] int32_t               rid,
                             [[maybe_unused]] ref_seqs_type       & ref_seqs,
                             [[maybe_unused]] int32_t               ref_start,
                             size_t                                 ref_length)
    {
        if (rid > -1 && ref_start > -1 &&       // read is mapped
            !cigar.empty() &&                   // alignment field was not empty
            !std::ranges::empty(get<1>(align))) // seq field was not empty
        {
            if constexpr (!detail::decays_to_ignore_v<ref_seqs_type>)
            {
                assert(static_cast<size_t>(ref_start + ref_length) <= std::ranges::size(ref_seqs[rid]));
                // copy over unaligned reference sequence part
                assign_unaligned(get<0>(align), ref_seqs[rid] | view::slice(ref_start, ref_start + ref_length));
            }
            else
            {
                using unaligned_t = remove_cvref_t<detail::unaligned_seq_t<decltype(get<0>(align))>>;
                auto dummy_seq    = view::repeat_n(value_type_t<unaligned_t>{}, ref_length)
                                  | std::view::transform(detail::access_restrictor_fn{});
                static_assert(std::Same<unaligned_t, decltype(dummy_seq)>,
                              "No reference information was given so the type of the first alignment tuple position"
                              "must have an unaligned sequence type of a dummy sequence ("
                              "view::repeat_n(dna5{}, size_t{}) | "
                              "std::view::transform(detail::access_restrictor_fn{}))");

                assign_unaligned(get<0>(align), dummy_seq); // assign dummy sequence
            }

            // insert gaps according to the cigar information
            detail::alignment_from_cigar(align, cigar);
        }
        else // not enough information for an alignment, assign an empty view/dummy_sequence
        {
            if constexpr (!detail::decays_to_ignore_v<ref_seqs_type>) // reference info given
            {
                assert(std::ranges::size(ref_seqs) > 0); // we assume that the given ref info is not empty
                assign_unaligned(get<0>(align), ref_seqs[0] | view::slice(0, 0));
            }
            else
            {
                using unaligned_t = remove_cvref_t<detail::unaligned_seq_t<decltype(get<0>(align))>>;
                assign_unaligned(get<0>(align), view::repeat_n(value_type_t<unaligned_t>{}, 0)
                                                | std::view::transform(detail::access_restrictor_fn{}));
            }
        }
    }

    template <typename stream_view_type>
    void read_field(stream_view_type && stream_view, detail::ignore_t const & SEQAN3_DOXYGEN_ONLY(target))
    {
        detail::consume(stream_view);
    }

    template <typename stream_view_type, std::ranges::ForwardRange target_range_type>
    void read_field(stream_view_type && stream_view, target_range_type & target)
    {
        if (!is_char<'*'>(*std::ranges::begin(stream_view)))
            std::ranges::copy(stream_view | view::char_to<value_type_t<target_range_type>>,
                              std::back_inserter(target));
        else
            std::ranges::next(std::ranges::begin(stream_view)); // skip '*'
    }

    template <typename stream_view_type, Arithmetic target_type>
    void read_field(stream_view_type && stream_view, target_type & target)
    {
        // unfortunately std::from_chars only accepts char const * so we need a buffer.
        auto [ignore, end] = std::ranges::copy(stream_view, buffer.data());
        (void) ignore;
        std::from_chars_result res = std::from_chars(buffer.begin(), end, target);

        if (res.ec == std::errc::invalid_argument || res.ptr != end)
            throw format_error{std::string("[CORRUPTED SAM FILE] The string '") + std::string(buffer.begin(), end) +
                                           "' could not be cast into type " +
                                           detail::get_display_name_v<target_type>.str()};

        if (res.ec == std::errc::result_out_of_range)
            throw format_error{std::string("[CORRUPTED SAM FILE] Casting '") + std::string(buffer.begin(), end) +
                                           "' into type " + detail::get_display_name_v<target_type>.str() +
                                           " would cause an overflow."};
    }

    template <typename stream_view_type, typename optional_value_type>
    void read_field(stream_view_type && stream_view, std::optional<optional_value_type> & target)
    {
        optional_value_type tmp;
        read_field(std::forward<stream_view_type>(stream_view), tmp);
        target = tmp;
    }

    template <typename stream_view_type, typename value_type>
    void read_sam_dict_vector(seqan3::detail::sam_tag_variant & variant,
                              stream_view_type && stream_view,
                              value_type value)
    {
        std::vector<value_type> tmp_vector;
        while (std::ranges::begin(stream_view) != ranges::end(stream_view)) // not fully consumed yet
        {
            read_field(stream_view | view::take_until(is_char<','>), value);
            tmp_vector.push_back(value);

            if (is_char<','>(*std::ranges::begin(stream_view)))
                std::ranges::next(std::ranges::begin(stream_view)); // skip ','
        }
        variant = std::move(tmp_vector);
    }

    template <typename stream_view_type>
    void read_field(stream_view_type && stream_view, sam_tag_dictionary & target)
    {
        /* Every SAM tag has the format "[TAG]:[TYPE_ID]:[VALUE]", where TAG is a two letter
           name tag which is converted to a unique integer identifier and TYPE_ID is one character in [A,i,Z,H,B,f]
           describing the type for the upcoming VALUES. If TYPE_ID=='B' it signals an array of comma separated
           VALUE's and the inner value type is identified by the character following ':', one of [cCsSiIf].
        */
        uint16_t tag = static_cast<uint16_t>(*std::ranges::begin(stream_view)) << 8;
        std::ranges::next(std::ranges::begin(stream_view)); // skip char read before
        tag += static_cast<uint16_t>(*std::ranges::begin(stream_view));
        std::ranges::next(std::ranges::begin(stream_view)); // skip char read before
        std::ranges::next(std::ranges::begin(stream_view)); // skip ':'
        char type_id = *std::ranges::begin(stream_view);
        std::ranges::next(std::ranges::begin(stream_view)); // skip char read before
        std::ranges::next(std::ranges::begin(stream_view)); // skip ':'

        switch (type_id)
        {
            case 'A' : // char
            {
                target[tag] = static_cast<char>(*std::ranges::begin(stream_view));
                std::ranges::next(std::ranges::begin(stream_view)); // skip char that has been read
                break;
            }
            case 'i' : // int32_t
            {
                int32_t tmp;
                read_field(stream_view, tmp);
                target[tag] = tmp;
                break;
            }
            case 'f' : // float
            {
                float tmp;
                read_field(stream_view, tmp);
                target[tag] = tmp;
                break;
            }
            case 'Z' : // string
            {
                target[tag] = std::string(stream_view);
                break;
            }
            case 'H' :
            {
                // TODO
                break;
            }
            case 'B' : // Array. Value type depends on second char [cCsSiIf]
            {
                char array_value_type_id = *std::ranges::begin(stream_view);
                std::ranges::next(std::ranges::begin(stream_view)); // skip char read before
                std::ranges::next(std::ranges::begin(stream_view)); // skip first ','

                switch (array_value_type_id)
                {
                    case 'c' : // int8_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, int8_t{});
                        break;
                    }
                    case 'C' : // uint8_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, uint8_t{});
                        break;
                    }
                    case 's' : // int16_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, int16_t{});
                        break;
                    }
                    case 'S' : // uint16_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, uint16_t{});
                        break;
                    }
                    case 'i' : // int32_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, int32_t{});
                        break;
                    }
                    case 'I' : // uint32_t
                    {
                        read_sam_dict_vector(target[tag], stream_view, uint32_t{});
                        break;
                    }
                    case 'f' : // float
                    {
                        read_sam_dict_vector(target[tag], stream_view, float{});
                        break;
                    }
                    default:
                        throw format_error{std::string("The first character in the numerical ") +
                                           "id of a SAM tag must be one of [cCsSiIf] but '" + array_value_type_id +
                                           "' was given."};
                }
                break;
            }
            default:
                throw format_error{std::string("The second character in the numerical id of a "
                                   "SAM tag must be one of [A,i,Z,H,B,f] but '") + type_id + "' was given."};
        }
    }

    template <typename stream_view_type, typename ref_ids_type, typename ref_seqs_type>
    void read_header(stream_view_type && stream_view,
                     alignment_file_header<ref_ids_type> & hdr,
                     ref_seqs_type & /*ref_id_to_pos_map*/)
    {
        auto parse_tag_value = [&stream_view, this] (auto & value) // helper function to parse the next tag value
        {
            detail::consume(stream_view | view::take_until_or_throw(is_char<':'>)); // skip tag name
            std::ranges::next(std::ranges::begin(stream_view));                     // skip ':'
            read_field(stream_view | view::take_until_or_throw(is_char<'\t'> || is_char<'\n'>), value);
        };

        // @HQ line
        // -------------------------------------------------------------------------------------------------------------
        parse_tag_value(hdr.format_version); // parse required VN (version) tag

        // The SO, SS and GO tag are optional and can appear in any order
        while (is_char<'\t'>(*std::ranges::begin(stream_view)))
        {
            std::ranges::next(std::ranges::begin(stream_view));              // skip tab
            std::string * who = std::addressof(hdr.grouping);

            if (is_char<'S'>(*std::ranges::begin(stream_view)))
            {
                std::ranges::next(std::ranges::begin(stream_view));          // skip S

                if (is_char<'O'>(*std::ranges::begin(stream_view)))          // SO (sorting) tag
                    who = std::addressof(hdr.sorting);
                else if (is_char<'S'>(*std::ranges::begin(stream_view)))     // SS (sub-order) tag
                    who = std::addressof(hdr.subsorting);
                else
                    throw format_error{std::string{"Illegal SAM header tag: S"} +
                                       std::string{static_cast<char>(*std::ranges::begin(stream_view))}};
            }
            else if (!is_char<'G'>(*std::ranges::begin(stream_view)))        // GO (grouping) tag
            {
                throw format_error{std::string{"Illegal SAM header tag in @HG starting with:"} +
                                   std::string{static_cast<char>(*std::ranges::begin(stream_view))}};
            }

            parse_tag_value(*who);
        }
        std::ranges::next(std::ranges::begin(stream_view));                  // skip newline

        // The rest of the header lines
        // -------------------------------------------------------------------------------------------------------------
        while (is_char<'@'>(*std::ranges::begin(stream_view)))
        {
            std::ranges::next(std::ranges::begin(stream_view));              // skip @

            if (is_char<'S'>(*std::ranges::begin(stream_view)))              // SQ (sequence dictionary) tag
            {
                ref_info_present_in_header = true;
                std::string id;
                std::tuple<int32_t, std::string> info{};

                parse_tag_value(id);                                         // parse required SN (sequence name) tag
                std::ranges::next(std::ranges::begin(stream_view));          // skip tab or newline
                parse_tag_value(get<0>(info));                               // parse required LN (length) tag

                if (is_char<'\t'>(*std::ranges::begin(stream_view)))         // read rest of the tags
                {
                    std::ranges::next(std::ranges::begin(stream_view));      // skip tab
                    read_field(stream_view | view::take_until_or_throw(is_char<'\n'>), get<1>(info));
                }
                std::ranges::next(std::ranges::begin(stream_view));          // skip newline

                /* If reference information were given, the ids exist and we can fill ref_dict directly.
                 * If not, wee need to update the ids first and fill the reference dictionary afterwards. */
                if constexpr (!detail::decays_to_ignore_v<ref_seqs_type>) // reference information given
                {
                    auto id_it = hdr.ref_dict.find(id);

                    if (id_it == hdr.ref_dict.end())
                        throw format_error{to_string("Unknown reference name '", id, "' found in SAM header ",
                                                     "(header.ref_ids(): ", hdr.ref_ids(), ").")};

                    auto & given_ref_info = hdr.ref_id_info[id_it->second];

                    if (std::get<0>(given_ref_info) != std::get<0>(info))
                        throw format_error{"Provided reference has unequal length as specified in the header."};

                    hdr.ref_id_info[id_it->second] = std::move(info);
                }
                else
                {
                    static_assert(!detail::is_type_specialisation_of_v<decltype(hdr.ref_ids()), std::deque>,
                                  "The range over reference ids must be of type std::deque such that "
                                  "pointers are not invalidated.");

                    hdr.ref_ids().push_back(id);
                    hdr.ref_id_info.push_back(info);
                    hdr.ref_dict[(hdr.ref_ids())[(hdr.ref_ids()).size() - 1]] = (hdr.ref_ids()).size() - 1;
                }
            }
            else if (is_char<'R'>(*std::ranges::begin(stream_view)))         // RG (read group) tag
            {
                std::pair<std::string, std::string> tmp{};

                parse_tag_value(get<0>(tmp));                                // read required ID tag

                if (is_char<'\t'>(*std::ranges::begin(stream_view)))         // read rest of the tags
                {
                    std::ranges::next(std::ranges::begin(stream_view));
                    read_field(stream_view | view::take_until_or_throw(is_char<'\n'>), get<1>(tmp));
                }
                std::ranges::next(std::ranges::begin(stream_view));          // skip newline

                hdr.read_groups.emplace_back(std::move(tmp));
            }
            else if (is_char<'P'>(*std::ranges::begin(stream_view)))         // PG (program) tag
            {
                typename alignment_file_header<ref_ids_type>::program_info_t tmp{};

                parse_tag_value(tmp.id);                                     // read required ID tag

                // The PN, CL, PP, DS, VN are optional tags and can be given in any order.
                while (is_char<'\t'>(*std::ranges::begin(stream_view)))
                {
                    std::ranges::next(std::ranges::begin(stream_view));      // skip tab
                    std::string * who = &tmp.version;

                    if (is_char<'P'>(*std::ranges::begin(stream_view)))
                    {
                        std::ranges::next(std::ranges::begin(stream_view));  // skip P

                        if (is_char<'N'>(*std::ranges::begin(stream_view)))  // PN (program name) tag
                            who = &tmp.name;
                        else                                                 // PP (previous program) tag
                            who = &tmp.previous;
                    }
                    else if (is_char<'C'>(*std::ranges::begin(stream_view))) // CL (command line) tag
                    {
                        who = &tmp.command_line_call;
                    }
                    else if (is_char<'D'>(*std::ranges::begin(stream_view))) // DS (description) tag
                    {
                        who = &tmp.description;
                    }
                    else if (!is_char<'V'>(*std::ranges::begin(stream_view))) // VN (version) tag
                    {
                        throw format_error{std::string{"Illegal SAM header tag starting with:"} +
                                           std::string{static_cast<char>(*std::ranges::begin(stream_view))}};
                    }

                    parse_tag_value(*who);
                }
                std::ranges::next(std::ranges::begin(stream_view));          // skip newline

                hdr.program_infos.emplace_back(std::move(tmp));
            }
            else if (is_char<'C'>(*std::ranges::begin(stream_view)))         // CO (comment) tag
            {
                std::string tmp;
                std::ranges::next(std::ranges::begin(stream_view)); // skip C
                std::ranges::next(std::ranges::begin(stream_view)); // skip O
                std::ranges::next(std::ranges::begin(stream_view)); // skip :
                read_field(stream_view | view::take_until_or_throw(is_char<'\n'>), tmp);
                std::ranges::next(std::ranges::begin(stream_view)); // skip newline

                hdr.comments.emplace_back(std::move(tmp));
            }
            else
            {
                throw format_error{std::string{"Illegal SAM header tag starting with:"} +
                                   std::string{static_cast<char>(*std::ranges::begin(stream_view))}};
            }
        }
    }
};

template <>
class alignment_file_output_format<format_sam>
{
public:
    using format_tag = format_sam;

    alignment_file_output_format()                                                 noexcept = default; 
    alignment_file_output_format(alignment_file_output_format const &)                      = delete;
    alignment_file_output_format & operator=(alignment_file_output_format const &)          = delete;
    alignment_file_output_format(alignment_file_output_format &&)                  noexcept = default; 
    alignment_file_output_format & operator=(alignment_file_output_format &&)      noexcept = default; 
    ~alignment_file_output_format()                                                noexcept = default; 

    template <typename stream_type,
              typename header_type,
              typename seq_type,
              typename id_type,
              typename ref_seq_type,
              typename ref_id_type,
              typename align_type,
              typename qual_type,
              typename mate_type,
              typename tag_dict_type>
    void write(stream_type                            &  stream,
               alignment_file_output_options const    &  options,
               header_type                            && header,
               seq_type                               && seq,
               qual_type                              && qual,
               id_type                                && id,
               int32_t                                   offset,
               ref_seq_type                           && SEQAN3_DOXYGEN_ONLY(ref_seq),
               ref_id_type                            && ref_id,
               std::optional<int32_t>                    ref_offset,
               align_type                             && align,
               uint16_t                                  flag,
               uint8_t                                   mapq,
               mate_type                              && mate,
               tag_dict_type                          && tag_dict,
               double                                    SEQAN3_DOXYGEN_ONLY(e_value),
               double                                    SEQAN3_DOXYGEN_ONLY(bit_score))
    {
        /* Note the following general things:
         *
         * - Given the SAM specifications, all fields may be empty
         *
         * - Arithmetic values default to 0 while all others default to '*'
         *
         * - Because of the former, arithmetic values can be directly streamed
         *   into 'stream' as operator<< is defined for all arithmetic types
         *   and the default value (0) is also the SAM default.
         *
         * - All other non-arithmetic values need to be checked for emptiness
         */

        // ---------------------------------------------------------------------
        // Type Requirements (as static asserts for user friendliness)
        // ---------------------------------------------------------------------
        static_assert((std::ranges::ForwardRange<seq_type>        &&
                      Alphabet<reference_t<seq_type>>),
                      "The seq object must be a std::ranges::ForwardRange over "
                      "letters that model seqan3::Alphabet.");

        static_assert((std::ranges::ForwardRange<id_type>         &&
                      Alphabet<reference_t<id_type>>),
                      "The id object must be a std::ranges::ForwardRange over "
                      "letters that model seqan3::Alphabet.");

        static_assert((std::ranges::ForwardRange<ref_seq_type>    &&
                      Alphabet<reference_t<ref_seq_type>>),
                      "The ref_seq object must be a std::ranges::ForwardRange "
                      "over letters that model seqan3::Alphabet.");

        if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
        {
            static_assert((std::ranges::ForwardRange<ref_id_type> ||
                           std::Integral<std::remove_reference_t<ref_id_type>> ||
                           detail::is_type_specialisation_of_v<remove_cvref_t<ref_id_type>, std::optional>),
                          "The ref_id object must be a std::ranges::ForwardRange "
                          "over letters that model seqan3::Alphabet.");

            if constexpr (std::Integral<remove_cvref_t<ref_id_type>> ||
                          detail::is_type_specialisation_of_v<remove_cvref_t<ref_id_type>, std::optional>)
                static_assert(!detail::decays_to_ignore_v<header_type>,
                              "If you give indices as reference id information the header must also be present.");
        }

        static_assert(TupleLike<remove_cvref_t<align_type>>,
                      "The align object must be a std::pair of two ranges whose "
                      "value_type is comparable to seqan3::gap");

        static_assert((std::tuple_size_v<remove_cvref_t<align_type>> == 2 &&
                       std::EqualityComparableWith<gap, reference_t<decltype(std::get<0>(align))>> &&
                       std::EqualityComparableWith<gap, reference_t<decltype(std::get<1>(align))>>),
                      "The align object must be a std::pair of two ranges whose "
                      "value_type is comparable to seqan3::gap");

        static_assert((std::ranges::ForwardRange<qual_type>       &&
                       Alphabet<reference_t<qual_type>>),
                      "The qual object must be a std::ranges::ForwardRange "
                      "over letters that model seqan3::Alphabet.");

        static_assert(TupleLike<remove_cvref_t<mate_type>>,
                      "The mate object must be a std::tuple of size 3 with "
                      "1) a std::ranges::ForwardRange with a value_type modelling seqan3::Alphabet, "
                      "2) a std::Integral or std::optional<std::Integral>, and "
                      "3) a std::Integral.");

        static_assert(((std::ranges::ForwardRange<decltype(std::get<0>(mate))>     ||
                        std::Integral<remove_cvref_t<decltype(std::get<0>(mate))>> ||
                        detail::is_type_specialisation_of_v<remove_cvref_t<decltype(std::get<0>(mate))>, std::optional>) &&
                      (std::Integral<remove_cvref_t<decltype(std::get<1>(mate))>> ||
                       detail::is_type_specialisation_of_v<remove_cvref_t<decltype(std::get<1>(mate))>, std::optional>) &&
                      std::Integral<remove_cvref_t<decltype(std::get<2>(mate))>>),
                      "The mate object must be a std::tuple of size 3 with "
                      "1) a std::ranges::ForwardRange with a value_type modelling seqan3::Alphabet, "
                      "2) a std::Integral or std::optional<std::Integral>, and "
                      "3) a std::Integral.");

        if constexpr (std::Integral<remove_cvref_t<decltype(std::get<0>(mate))>> ||
                      detail::is_type_specialisation_of_v<remove_cvref_t<decltype(std::get<0>(mate))>, std::optional>)
            static_assert(!detail::decays_to_ignore_v<header_type>,
                          "If you give indices as mate reference id information the header must also be present.");

        static_assert(std::Same<remove_cvref_t<tag_dict_type>, sam_tag_dictionary>,
                      "The tag_dict object must be of type seqan3::sam_tag_dictionary.");

        // ---------------------------------------------------------------------
        // logical Requirements
        // ---------------------------------------------------------------------
        if constexpr (!detail::decays_to_ignore_v<header_type> &&
                      !detail::decays_to_ignore_v<ref_id_type> &&
                      !std::Integral<std::remove_reference_t<ref_id_type>> &&
                      !detail::is_type_specialisation_of_v<std::remove_reference_t<ref_id_type>, std::optional>)
        {
            static_assert(ImplicitlyConvertibleTo<ref_id_type &, typename decltype(header.ref_dict)::key_type>,
                      "The ref_id type is not convertible to the reference id information stored in the "
                      "reference dictionary of the header object.");

            if (options.sam_require_header && !std::ranges::empty(ref_id))
            {
                if ((header.ref_dict).count(ref_id) == 0) // no reference id matched
                    throw format_error{std::string("The ref_id '") + std::string(ref_id) +
                                       "' was not in the list of references"};
            }
        }

        if (ref_offset.has_value() && (ref_offset.value() + 1) < 0)
            throw format_error{"The ref_offset object must be an std::Integral >= 0."};

        // ---------------------------------------------------------------------
        // Writing the Header on first call
        // ---------------------------------------------------------------------
        if constexpr (!detail::decays_to_ignore_v<header_type>)
        {
            if (options.sam_require_header && !written_header)
            {
                write_header(stream, options, header);
                written_header = true;
            }
        }

        // ---------------------------------------------------------------------
        // Writing the Record
        // ---------------------------------------------------------------------
        seqan3::ostreambuf_iterator stream_it{stream};
        char const separator{'\t'};

        write_range(stream_it, std::forward<id_type>(id));

        stream << separator;

        stream << flag << separator;

        if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
        {
            if constexpr (std::Integral<std::remove_reference_t<ref_id_type>>)
            {
                write_range(stream_it, (header.ref_ids())[ref_id]);
            }
            else if constexpr (detail::is_type_specialisation_of_v<std::remove_reference_t<ref_id_type>, std::optional>)
            {
                if (ref_id.has_value())
                    write_range(stream_it, (header.ref_ids())[ref_id.value()]);
                else
                    stream << '*';
            }
            else
            {
                write_range(stream_it, std::forward<ref_id_type>(ref_id));
            }
        }
        else
        {
            stream << '*';
        }

        stream << separator;

        // SAM is 1 based, 0 indicates unmapped read if optional is not set
        stream << (ref_offset.value_or(-1) + 1) << separator;

        stream << static_cast<unsigned>(mapq) << separator;

        if (!std::ranges::empty(get<0>(align)) && !std::ranges::empty(get<1>(align)))
        {
            // compute possible distance from alignment end to sequence end
            // which indicates soft clipping at the end.
            // This should be replace by a free count_gaps function for
            // aligned sequences which is more efficient if possible.
            size_t off_end{std::ranges::size(seq) - offset};
            for (auto chr : get<1>(align))
                if (chr == gap{})
                    ++off_end;
            off_end -= std::ranges::size(get<1>(align));

            write_range(stream_it, detail::get_cigar_string(std::forward<align_type>(align), offset, off_end));
        }
        else
        {
            stream << '*';
        }

        stream << separator;

        if constexpr (std::Integral<std::remove_reference_t<decltype(get<0>(mate))>>)
        {
            write_range(stream_it, (header.ref_ids())[get<0>(mate)]);
        }
        else if constexpr (detail::is_type_specialisation_of_v<std::remove_reference_t<decltype(get<0>(mate))>, std::optional>)
        {
            if (get<0>(mate).has_value())
                // value_or(0) instead of value() (which is equivalent here) as a
                // workaround for a ubsan false-positive in GCC8: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90058
                write_range(stream_it, header.ref_ids()[get<0>(mate).value_or(0)]);
            else
                stream << '*';
        }
        else
        {
            write_range(stream_it, get<0>(mate));
        }

        stream << separator;

        if constexpr (detail::is_type_specialisation_of_v<remove_cvref_t<decltype(get<1>(mate))>, std::optional>)
        {
            // SAM is 1 based, 0 indicates unmapped read if optional is not set
            stream << (get<1>(mate).value_or(-1) + 1) << separator;
        }
        else
        {
            stream << get<1>(mate) << separator;
        }

        stream << get<2>(mate) << separator;

        write_range(stream_it, std::forward<seq_type>(seq));

        stream << separator;

        write_range(stream_it, std::forward<qual_type>(qual));

        write_tag_fields(stream, std::forward<tag_dict_type>(tag_dict), separator);

        detail::write_eol(stream_it, options.add_carriage_return);
    }

protected:
    static constexpr char format_version[4] = "1.6";
    bool written_header{false};

    template <typename stream_it_t, typename field_type>
        requires std::ranges::ForwardRange<field_type>
    void write_range(stream_it_t & stream_it, field_type && field_value)
    {
        if (std::ranges::empty(field_value))
            stream_it = '*';
        else
            std::ranges::copy(field_value | view::to_char, stream_it);
    }

    template <typename stream_it_t>
    void write_range(stream_it_t & stream_it, char const * const field_value)
    {
        write_range(stream_it, std::string_view{field_value});
    }

    template <typename stream_t, Arithmetic field_type>
    void write_field(stream_t & stream, field_type field_value)
    {
        // TODO: replace this with to_chars for efficiency
        if constexpr (std::Same<field_type, int8_t> || std::Same<field_type, uint8_t>)
            stream << static_cast<int16_t>(field_value);
        else
            stream << field_value;
    }

    template <typename stream_t>
    void write_tag_fields(stream_t & stream, sam_tag_dictionary const & tag_dict, char const separator)
    {
        auto stream_variant_fn = [this, &stream] (auto && arg) // helper to print an std::variant
        {
            using T = remove_cvref_t<decltype(arg)>;

            if constexpr (!Container<T> || std::Same<T, std::string>)
            {
                stream << arg;
            }
            else
            {
                if (arg.begin() != arg.end())
                {
                    for (auto it = arg.begin(); it != (arg.end() - 1); ++it)
                    {
                        write_field(stream, *it);
                        stream << ',';
                    }

                    write_field(stream, *(arg.end() - 1)); // write last value without trailing ','
                }
            }
        };

        for (auto & [tag, variant] : tag_dict)
        {
            stream << separator;

            char char0 = tag / 256;
            char char1 = tag % 256;

            stream << char0 << char1 << ':' << detail::sam_tag_type_char[variant.index()] << ':';

            if (detail::sam_tag_type_char_extra[variant.index()] != '\0')
                stream << detail::sam_tag_type_char_extra[variant.index()] << ',';

            std::visit(stream_variant_fn, variant);
        }
    }

    template <typename stream_t, typename ref_ids_type>
    void write_header(stream_t & stream,
                      alignment_file_output_options const & options,
                      alignment_file_header<ref_ids_type> & header)
    {
        // -----------------------------------------------------------------
        // Check Header
        // -----------------------------------------------------------------

        // (@HD) Check header line
        // The format version string will be taken from the local member variable
        if (!header.sorting.empty() &&
            !(header.sorting == "unknown"   ||
              header.sorting == "unsorted"  ||
              header.sorting == "queryname" ||
              header.sorting == "coordinate" ))
            throw format_error{"SAM format error: The header.sorting member must be "
                               "one of [unknown, unsorted, queryname, coordinate]."};

        if (!header.grouping.empty() &&
            !(header.grouping == "none"   ||
              header.grouping == "query"  ||
              header.grouping == "reference"))
            throw format_error{"SAM format error: The header.grouping member must be "
                               "one of [none, query, reference]."};

        // (@SQ) Check Reference Sequence Dictionary lines

        // TODO

        // - sorting order be one of ...
        // - grouping can be one of ...
        // - reference names must be unique
        // - ids of read groups must be unique
        // - program ids need to be unique
        // many more small semantic things, like fits REGEX

        // -----------------------------------------------------------------
        // Write Header
        // -----------------------------------------------------------------
        seqan3::ostreambuf_iterator stream_it{stream};

        // (@HD) Write header line [required].
        stream << "@HD\tVN:";
        stream << format_sam::format_version;

        if (!header.sorting.empty())
            stream << "\tSO:" << header.sorting;

        if (!header.subsorting.empty())
            stream << "\tSS:" << header.subsorting;

        if (!header.grouping.empty())
            stream << "\tGO:" << header.grouping;

        detail::write_eol(stream_it, options.add_carriage_return);

        // (@SQ) Write Reference Sequence Dictionary lines [required].
        for (auto const & [ref_name, ref_info] : std::view::zip(header.ref_ids(), header.ref_id_info))
        {
            stream << "@SQ\tSN:";

            std::ranges::copy(ref_name, stream_it);

            stream << "\tLN:" << get<0>(ref_info);

            if (!get<1>(ref_info).empty())
                stream << "\t" << get<1>(ref_info);

            detail::write_eol(stream_it, options.add_carriage_return);
        }

        // Write read group (@RG) lines if specified.
        for (auto const & read_group : header.read_groups)
        {
            stream << "@RG"
                   << "\tID:" << get<0>(read_group);

            if (!get<1>(read_group).empty())
                stream << "\t" << get<1>(read_group);

            detail::write_eol(stream_it, options.add_carriage_return);
        }

        // Write program (@PG) lines if specified.
        for (auto const & program : header.program_infos)
        {
            stream << "@PG"
                   << "\tID:" << program.id;

            if (!program.name.empty())
                stream << "\tPN:" << program.name;

            if (!program.command_line_call.empty())
                stream << "\tCL:" << program.command_line_call;

            if (!program.previous.empty())
                stream << "\tPP:" << program.previous;

            if (!program.description.empty())
                stream << "\tDS:" << program.description;

            if (!program.version.empty())
                stream << "\tVN:" << program.version;

            detail::write_eol(stream_it, options.add_carriage_return);
        }

        // Write comment (@CO) lines if specified.
        for (auto const & comment : header.comments)
        {
            stream << "@CO\t" << comment;
            detail::write_eol(stream_it, options.add_carriage_return);
        }
    }
};

} // namespace seqan3::detail