format_fasta.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <algorithm>
#include <iterator>
#include <ranges>
#include <string>
#include <string_view>
#include <vector>
 
#include <seqan3/alphabet/adaptation/char.hpp>
#include <seqan3/alphabet/concept.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/alphabet/quality/aliases.hpp>
#include <seqan3/alphabet/views/char_to.hpp>
#include <seqan3/alphabet/views/to_char.hpp>
#include <seqan3/core/range/detail/misc.hpp>
#include <seqan3/core/range/type_traits.hpp>
#include <seqan3/io/detail/ignore_output_iterator.hpp>
#include <seqan3/io/detail/misc.hpp>
#include <seqan3/io/sequence_file/input_format_concept.hpp>
#include <seqan3/io/sequence_file/input_options.hpp>
#include <seqan3/io/sequence_file/output_format_concept.hpp>
#include <seqan3/io/sequence_file/output_options.hpp>
#include <seqan3/io/stream/detail/fast_ostreambuf_iterator.hpp>
#include <seqan3/io/views/detail/istreambuf_view.hpp>
#include <seqan3/io/views/detail/take_line_view.hpp>
#include <seqan3/io/views/detail/take_until_view.hpp>
#include <seqan3/utility/char_operations/predicate.hpp>
#include <seqan3/utility/detail/type_name_as_string.hpp>
 
namespace seqan3
{
 
class format_fasta
{
public:
    format_fasta() noexcept = default;                                 
    format_fasta(format_fasta const &) noexcept = default;             
    format_fasta & operator=(format_fasta const &) noexcept = default; 
    format_fasta(format_fasta &&) noexcept = default;                  
    format_fasta & operator=(format_fasta &&) noexcept = default;      
    ~format_fasta() noexcept = default;                                
 
 
    static inline std::vector<std::string> file_extensions{
        {"fasta"},
        {"fa"},
        {"fna"},
        {"ffn"},
        {"faa"},
        {"frn"},
        {"fas"},
    };
 
protected:
    template <typename stream_type, // constraints checked by file
              typename legal_alph_type,
              typename stream_pos_type,
              typename seq_type, // other constraints checked inside function
              typename id_type,
              typename qual_type>
    void read_sequence_record(stream_type & stream,
                              sequence_file_input_options<legal_alph_type> const & options,
                              stream_pos_type & position_buffer,
                              seq_type & sequence,
                              id_type & id,
                              qual_type & SEQAN3_DOXYGEN_ONLY(qualities))
    {
        // Store current position in buffer
        // Must happen before constructing the view.
        // With libc++, tellg invalidates the I/O buffer.
        position_buffer = stream.tellg();
 
        auto stream_view = detail::istreambuf(stream);
 
        // ID
        read_id(stream_view, options, id);
 
        // Sequence
        read_seq(stream_view, options, sequence);
    }
 
    template <typename stream_type, // constraints checked by file
              typename seq_type,    // other constraints checked inside function
              typename id_type,
              typename qual_type>
    void write_sequence_record(stream_type & stream,
                               sequence_file_output_options const & options,
                               seq_type && sequence,
                               id_type && id,
                               qual_type && SEQAN3_DOXYGEN_ONLY(qualities))
    {
        seqan3::detail::fast_ostreambuf_iterator stream_it{*stream.rdbuf()};
 
        // ID
        if constexpr (detail::decays_to_ignore_v<id_type>)
        {
            throw std::logic_error{"The ID field may not be set to ignore when writing FASTA files."};
        }
        else
        {
            if (std::ranges::empty(id)) //[[unlikely]]
                throw std::runtime_error{"The ID field may not be empty when writing FASTA files."};
 
            write_id(stream_it, options, id);
        }
 
        // Sequence
        if constexpr (detail::decays_to_ignore_v<seq_type>) // sequence
        {
            throw std::logic_error{
                "The SEQ and SEQ_QUAL fields may not both be set to ignore when writing FASTA files."};
        }
        else
        {
            if (std::ranges::empty(sequence)) //[[unlikely]]
                throw std::runtime_error{"The SEQ field may not be empty when writing FASTA files."};
 
            write_seq(stream_it, options, sequence);
        }
    }
 
private:
    template <typename stream_view_t, typename seq_legal_alph_type, typename id_type>
    void
    read_id(stream_view_t & stream_view, sequence_file_input_options<seq_legal_alph_type> const & options, id_type & id)
    {
        auto const is_id = is_char<'>'> || is_char<';'>;
 
        if (!is_id(*begin(stream_view)))
            throw parse_error{std::string{"Expected to be on beginning of ID, but "} + is_id.msg
                              + " evaluated to false on " + detail::make_printable(*begin(stream_view))};
 
        if constexpr (detail::decays_to_ignore_v<id_type>) // Skip the ID, it is not requested by the user.
        {
            detail::consume(stream_view | detail::take_line_or_throw);
        }
        else // read ID
        {
            if (options.truncate_ids)
            {
#if SEQAN3_WORKAROUND_VIEW_PERFORMANCE
                auto it = stream_view.begin();
                auto e = stream_view.end();
                ++it; // already checked `is_id`
 
                if (options.fasta_ignore_blanks_before_id)
                {
                    for (; (it != e) && (is_blank)(*it); ++it) // skip leading ' '
                    {}
                }
 
                bool at_delimiter = false;
                for (; it != e; ++it)
                {
                    if ((is_cntrl || is_blank)(*it))
                    {
                        at_delimiter = true;
                        break;
                    }
                    id.push_back(assign_char_to(*it, std::ranges::range_value_t<id_type>{}));
                }
 
                if (!at_delimiter)
                    throw unexpected_end_of_input{"FASTA ID line did not end in newline."};
 
                for (; (it != e) && ((!is_char<'\n'>)(*it)); ++it)
                {}
 
#else  // ↑↑↑ WORKAROUND | ORIGINAL ↓↓↓
                if (options.fasta_ignore_blanks_before_id)
                {
                    std::ranges::copy(
                        stream_view | std::views::drop(1)                       // skip leading '>' or ';'
                            | std::views::drop_while(is_blank)                  // skip leading ' '
                            | detail::take_until_or_throw(is_cntrl || is_blank) // read ID until delimiter…
                            | views::char_to<std::ranges::range_value_t<id_type>>,
                        std::back_inserter(id)); // … ^A is old delimiter
                }
                else
                {
                    std::ranges::copy(
                        stream_view | std::views::drop(1)                       // skip leading '>' or ';'
                            | detail::take_until_or_throw(is_cntrl || is_blank) // read ID until delimiter…
                            | views::char_to<std::ranges::range_value_t<id_type>>,
                        std::back_inserter(id)); // … ^A is old delimiter
                }
 
                // consume rest of line
                detail::consume(stream_view | detail::take_line_or_throw);
#endif // SEQAN3_WORKAROUND_VIEW_PERFORMANCE
            }
            else // options.truncate_ids
            {
#if SEQAN3_WORKAROUND_VIEW_PERFORMANCE
                auto it = stream_view.begin();
                auto e = stream_view.end();
                ++it; // skip leading '>' or ';'
 
                if (options.fasta_ignore_blanks_before_id)
                {
                    for (; (it != e) && (is_blank)(*it); ++it) // skip leading ' '
                    {}
                }
 
                bool at_delimiter = false;
                for (; it != e; ++it)
                {
                    if ((is_char<'\n'>)(*it))
                    {
                        at_delimiter = true;
                        break;
                    }
                    id.push_back(assign_char_to(*it, std::ranges::range_value_t<id_type>{}));
                }
 
                if (!at_delimiter)
                    throw unexpected_end_of_input{"FASTA ID line did not end in newline."};
 
#else  // ↑↑↑ WORKAROUND | ORIGINAL ↓↓↓
                if (options.fasta_ignore_blanks_before_id)
                {
                    std::ranges::copy(stream_view | detail::take_line_or_throw // read line
                                          | std::views::drop(1)                // skip leading '>' or ';'
                                          | std::views::drop_while(is_blank)   // skip leading ' '
                                          | views::char_to<std::ranges::range_value_t<id_type>>,
                                      std::back_inserter(id));
                }
                else
                {
                    std::ranges::copy(stream_view | detail::take_line_or_throw // read line
                                          | std::views::drop(1)                // skip leading '>' or ';'
                                          | views::char_to<std::ranges::range_value_t<id_type>>,
                                      std::back_inserter(id));
                }
#endif // SEQAN3_WORKAROUND_VIEW_PERFORMANCE
            }
        }
    }
 
    template <typename stream_view_t, typename seq_legal_alph_type, typename seq_type>
    void read_seq(stream_view_t & stream_view, sequence_file_input_options<seq_legal_alph_type> const &, seq_type & seq)
    {
        constexpr auto is_id = is_char<'>'> || is_char<';'>;
 
        if constexpr (!detail::decays_to_ignore_v<seq_type>)
        {
            constexpr auto is_legal_alph = char_is_valid_for<seq_legal_alph_type>;
 
#if SEQAN3_WORKAROUND_VIEW_PERFORMANCE
            auto it = stream_view.begin();
            auto e = stream_view.end();
 
            if (it == e)
                throw unexpected_end_of_input{"No sequence information given!"};
 
            for (; (it != e) && ((!is_id)(*it)); ++it)
            {
                if ((is_space || is_digit)(*it))
                {
                    continue;
                }
                else if (is_legal_alph(*it))
                {
                    seq.push_back(assign_char_to(*it, std::ranges::range_value_t<seq_type>{}));
                }
                else
                {
 
                    throw parse_error{std::string{"Encountered an unexpected letter: "} + "char_is_valid_for<"
                                      + detail::type_name_as_string<seq_legal_alph_type>
                                      + "> evaluated to false on " + detail::make_printable(*it)};
                }
            }
 
#else  // ↑↑↑ WORKAROUND | ORIGINAL ↓↓↓
 
            if (std::ranges::begin(stream_view) == std::ranges::end(stream_view))
                throw unexpected_end_of_input{"No sequence information given!"};
 
            std::ranges::copy(
                stream_view | detail::take_until(is_id)           // until next header (or end)
                    | std::views::filter(!(is_space || is_digit)) // ignore whitespace and numbers
                    | std::views::transform(
                        [is_legal_alph](char const c)
                        {
                            if (!is_legal_alph(c))
                            {
                                throw parse_error{std::string{"Encountered an unexpected letter: "}
                                                  + "char_is_valid_for<"
                                                  + detail::type_name_as_string<seq_legal_alph_type>
                                                  + "> evaluated to false on " + detail::make_printable(c)};
                            }
                            return c;
                        })                                                  // enforce legal alphabet
                    | views::char_to<std::ranges::range_value_t<seq_type>>, // convert to actual target alphabet
                std::back_inserter(seq));
#endif // SEQAN3_WORKAROUND_VIEW_PERFORMANCE
        }
        else
        {
            detail::consume(stream_view | detail::take_until(is_id));
        }
    }
 
    template <typename stream_it_t, typename id_type>
    void write_id(stream_it_t & stream_it, sequence_file_output_options const & options, id_type && id)
    {
        if (options.fasta_legacy_id_marker)
            stream_it = ';';
        else
            stream_it = '>';
 
        if (options.fasta_blank_before_id)
            stream_it = ' ';
 
        stream_it.write_range(id);
        stream_it.write_end_of_line(options.add_carriage_return);
    }
 
    template <typename stream_it_t, typename seq_type>
    void write_seq(stream_it_t & stream_it, sequence_file_output_options const & options, seq_type && seq)
    {
        auto char_sequence = seq | views::to_char;
 
        if (options.fasta_letters_per_line > 0)
        {
            /* Using `views::interleave` is probably the way to go but that needs performance-tuning.*/
            auto it = std::ranges::begin(char_sequence);
            auto end = std::ranges::end(char_sequence);
 
            while (it != end)
            {
                /* Note: This solution is slightly suboptimal for sized but non-random-access ranges.*/
                auto current_end = it;
                size_t steps = std::ranges::advance(current_end, options.fasta_letters_per_line, end);
                using subrange_t = std::ranges::subrange<decltype(it), decltype(it), std::ranges::subrange_kind::sized>;
                it = stream_it.write_range(subrange_t{it, current_end, (options.fasta_letters_per_line - steps)});
                stream_it.write_end_of_line(options.add_carriage_return);
            }
        }
        else
        {
            stream_it.write_range(char_sequence);
            stream_it.write_end_of_line(options.add_carriage_return);
        }
    }
};
 
} // namespace seqan3