31namespace seqan3::detail
48 format_sam_base() =
default;
49 format_sam_base(format_sam_base
const &) =
default;
50 format_sam_base & operator=(format_sam_base
const &) =
default;
51 format_sam_base(format_sam_base &&) =
default;
52 format_sam_base & operator=(format_sam_base &&) =
default;
53 ~format_sam_base() =
default;
58 static constexpr std::array format_version{
'1',
'.',
'6'};
64 bool header_was_written{
false};
67 bool ref_info_present_in_header{
false};
69 template <
typename ref_
id_type,
typename ref_
id_tmp_type,
typename header_type,
typename ref_seqs_type>
70 void check_and_assign_ref_id(ref_id_type &
ref_id,
71 ref_id_tmp_type & ref_id_tmp,
77 template <
typename stream_view_type, std::ranges::forward_range target_range_type>
78 void read_forward_range_field(stream_view_type && stream_view, target_range_type & target);
80 template <std::ranges::forward_range target_range_type>
81 void read_forward_range_field(
std::string_view const str, target_range_type & target);
83 template <arithmetic arithmetic_target_type>
84 void read_arithmetic_field(
std::string_view const & str, arithmetic_target_type & arithmetic_target);
86 template <
typename stream_view_type,
typename ref_
ids_type,
typename ref_seqs_type,
typename seq_legal_alph_type>
87 void read_header(stream_view_type && stream_view,
88 sam_file_header<ref_ids_type> & hdr,
90 sam_file_input_options<seq_legal_alph_type>
const & options);
92 template <
typename stream_t,
typename header_type>
93 void write_header(stream_t & stream, sam_file_output_options
const & options, header_type & header);
106template <
typename ref_
id_type,
typename ref_
id_tmp_type,
typename header_type,
typename ref_seqs_type>
107inline void format_sam_base::check_and_assign_ref_id(ref_id_type &
ref_id,
108 ref_id_tmp_type & ref_id_tmp,
109 header_type & header,
112 if (!std::ranges::empty(ref_id_tmp))
114 auto search = header.ref_dict.find(ref_id_tmp);
116 if (search == header.ref_dict.end())
118 if constexpr (detail::decays_to_ignore_v<ref_seqs_type>)
120 if (ref_info_present_in_header)
122 throw format_error{
"Unknown reference id found in record which is not present in the header."};
126 header.ref_ids().push_back(ref_id_tmp);
127 auto pos = std::ranges::size(header.ref_ids()) - 1;
128 header.ref_dict[header.ref_ids()[pos]] = pos;
134 throw format_error{
"Unknown reference id found in record which is not present in the given ids."};
147inline int32_t format_sam_base::soft_clipping_at_front(
std::vector<cigar> const & cigar_vector)
const
152 auto soft_clipping_at = [&](
size_t const index)
154 return cigar_vector[index] ==
'S'_cigar_operation;
157 auto hard_clipping_at = [&](
size_t const index)
159 return cigar_vector[index] ==
'H'_cigar_operation;
162 auto vector_size_at_least = [&](
size_t const min_size)
164 return cigar_vector.
size() >= min_size;
167 auto cigar_count_at = [&](
size_t const index)
169 return get<0>(cigar_vector[index]);
173 if (vector_size_at_least(1) && soft_clipping_at(0))
174 sc_front = cigar_count_at(0);
175 else if (vector_size_at_least(2) && hard_clipping_at(0) && soft_clipping_at(1))
176 sc_front = cigar_count_at(1);
188template <
typename stream_view_type, std::ranges::forward_range target_range_type>
189inline void format_sam_base::read_forward_range_field(stream_view_type && stream_view, target_range_type & target)
191 using target_range_value_t = std::ranges::range_value_t<target_range_type>;
192 using begin_iterator_t = std::ranges::iterator_t<stream_view_type>;
193 using end_iterator_t = std::ranges::sentinel_t<stream_view_type>;
200 if (
char c = *it; !(++it == std::ranges::end(stream_view) && c ==
'*'))
203 std::ranges::copy(std::ranges::subrange<begin_iterator_t, end_iterator_t>{it, std::ranges::end(stream_view)}
204 | views::char_to<target_range_value_t>,
216template <std::ranges::forward_range target_range_type>
217inline void format_sam_base::read_forward_range_field(
std::string_view const str, target_range_type & target)
219 if (str.
size() == 1 && str[0] ==
'*')
222 if constexpr (std::assignable_from<target_range_type, std::string_view>)
228 target.resize(str.
size());
229 for (
size_t i = 0; i < str.
size(); ++i)
230 target[i] =
assign_char_to(str[i], std::ranges::range_value_t<target_range_type>{});
243template <arithmetic arithmetic_target_type>
244inline void format_sam_base::read_arithmetic_field(
std::string_view const & str,
245 arithmetic_target_type & arithmetic_target)
249 if (res.ec == std::errc::invalid_argument || res.ptr != str.
end())
251 +
"' could not be cast into type " + detail::type_name_as_string<arithmetic_target_type>};
253 if (res.ec == std::errc::result_out_of_range)
255 +
"' into type " + detail::type_name_as_string<arithmetic_target_type>
256 +
" would cause an overflow."};
278template <
typename stream_view_type,
typename ref_
ids_type,
typename ref_seqs_type,
typename seq_legal_alph_type>
279inline void format_sam_base::read_header(stream_view_type && stream_view,
280 sam_file_header<ref_ids_type> & hdr,
282 sam_file_input_options<seq_legal_alph_type>
const & options)
285 auto end = std::ranges::end(stream_view);
288 auto make_tag = [](uint8_t char1, uint8_t char2)
constexpr
290 return static_cast<uint16_t
>(char1) | (
static_cast<uint16_t
>(char2) << CHAR_BIT);
295 auto parse_and_make_tag = [&]()
301 return make_tag(raw_tag[0], raw_tag[1]);
304 auto take_until_predicate = [&it, &string_buffer](
auto const & predicate)
306 string_buffer.clear();
307 while (!predicate(*it))
309 string_buffer.push_back(*it);
314 auto skip_until_predicate = [&it](
auto const & predicate)
316 while (!predicate(*it))
320 auto copy_next_tag_value_into_buffer = [&]()
322 skip_until_predicate(
is_char<
':'>);
341 read_forward_range_field(string_buffer, value);
344 auto print_cerr_of_unspported_tag = [&options](
char const *
const header_tag,
std::array<char, 2> raw_tag)
346 if (options.stream_warnings_to ==
nullptr)
349 *options.stream_warnings_to <<
"Unsupported SAM header tag in @" << header_tag <<
": " << raw_tag[0]
350 << raw_tag[1] <<
'\n';
353 while (it != end &&
is_char<
'@'>(*it))
357 switch (parse_and_make_tag())
359 case make_tag(
'H',
'D'):
367 switch (parse_and_make_tag())
369 case make_tag(
'V',
'N'):
374 case make_tag(
'S',
'O'):
379 case make_tag(
'S',
'S'):
384 case make_tag(
'G',
'O'):
391 print_cerr_of_unspported_tag(
"HD", raw_tag);
395 if (header_entry !=
nullptr)
397 copy_next_tag_value_into_buffer();
398 read_forward_range_field(string_buffer, *header_entry);
405 if (hdr.format_version.empty())
406 throw format_error{
std::string{
"The required VN tag in @HD is missing."}};
411 case make_tag(
'S',
'Q'):
413 ref_info_present_in_header =
true;
414 std::ranges::range_value_t<
decltype(hdr.ref_ids())>
id;
423 switch (parse_and_make_tag())
425 case make_tag(
'S',
'N'):
427 copy_next_tag_value_into_buffer();
428 read_forward_range_field(string_buffer,
id);
431 case make_tag(
'L',
'N'):
433 int32_t sequence_length_tmp{};
434 copy_next_tag_value_into_buffer();
435 read_arithmetic_field(string_buffer, sequence_length_tmp);
436 sequence_length = sequence_length_tmp;
441 parse_and_append_unhandled_tag_to_string(get<1>(info), raw_tag);
448 throw format_error{
std::string{
"The required SN tag in @SQ is missing."}};
449 if (!sequence_length.has_value())
450 throw format_error{
std::string{
"The required LN tag in @SQ is missing."}};
451 if (sequence_length.value() <= 0)
452 throw format_error{
std::string{
"The value of LN in @SQ must be positive."}};
454 get<0>(info) = sequence_length.value();
457 if constexpr (!detail::decays_to_ignore_v<ref_seqs_type>)
459 auto id_it = hdr.ref_dict.
find(
id);
461 if (id_it == hdr.ref_dict.end())
462 throw format_error{detail::to_string(
"Unknown reference name '",
464 "' found in SAM header ",
465 "(header.ref_ids(): ",
469 auto & given_ref_info = hdr.ref_id_info[id_it->second];
471 if (std::get<0>(given_ref_info) != std::get<0>(info))
472 throw format_error{
"Provided and header-based reference length differ."};
474 hdr.ref_id_info[id_it->second] = std::move(info);
478 static_assert(!detail::is_type_specialisation_of_v<
decltype(hdr.ref_ids()),
std::deque>,
479 "The range over reference ids must be of type std::deque such that pointers are not "
482 hdr.ref_ids().push_back(
id);
483 hdr.ref_id_info.push_back(info);
484 hdr.ref_dict[(hdr.ref_ids())[(hdr.ref_ids()).
size() - 1]] = (hdr.ref_ids()).
size() - 1;
489 case make_tag(
'R',
'G'):
498 switch (parse_and_make_tag())
500 case make_tag(
'I',
'D'):
502 copy_next_tag_value_into_buffer();
503 read_forward_range_field(string_buffer, get<0>(tmp));
508 parse_and_append_unhandled_tag_to_string(get<1>(tmp), raw_tag);
514 if (get<0>(tmp).
empty())
515 throw format_error{
std::string{
"The required ID tag in @RG is missing."}};
517 hdr.read_groups.emplace_back(std::move(tmp));
521 case make_tag(
'P',
'G'):
531 switch (parse_and_make_tag())
533 case make_tag(
'I',
'D'):
538 case make_tag(
'P',
'N'):
543 case make_tag(
'P',
'P'):
548 case make_tag(
'C',
'L'):
553 case make_tag(
'D',
'S'):
558 case make_tag(
'V',
'N'):
565 print_cerr_of_unspported_tag(
"PG", raw_tag);
569 if (program_info_entry !=
nullptr)
571 copy_next_tag_value_into_buffer();
572 read_forward_range_field(string_buffer, *program_info_entry);
580 throw format_error{
std::string{
"The required ID tag in @PG is missing."}};
582 hdr.program_infos.emplace_back(std::move(tmp));
586 case make_tag(
'C',
'O'):
590 take_until_predicate(
is_char<
'\n'>);
591 read_forward_range_field(string_buffer, tmp);
593 hdr.comments.emplace_back(std::move(tmp));
598 throw format_error{
std::string{
"Illegal SAM header tag starting with:"} + *it};
619template <
typename stream_t,
typename header_type>
621format_sam_base::write_header(stream_t & stream, sam_file_output_options
const & options, header_type & header)
623 if constexpr (!detail::decays_to_ignore_v<header_type>)
631 if (!header.sorting.empty()
632 && !(header.sorting ==
"unknown" || header.sorting ==
"unsorted" || header.sorting ==
"queryname"
633 || header.sorting ==
"coordinate"))
634 throw format_error{
"SAM format error: The header.sorting member must be "
635 "one of [unknown, unsorted, queryname, coordinate]."};
637 if (!header.grouping.empty()
638 && !(header.grouping ==
"none" || header.grouping ==
"query" || header.grouping ==
"reference"))
639 throw format_error{
"SAM format error: The header.grouping member must be "
640 "one of [none, query, reference]."};
659 stream <<
"@HD\tVN:";
662 if (!header.sorting.empty())
663 stream <<
"\tSO:" << header.sorting;
665 if (!header.subsorting.empty())
666 stream <<
"\tSS:" << header.subsorting;
668 if (!header.grouping.empty())
669 stream <<
"\tGO:" << header.grouping;
671 detail::write_eol(stream_it, options.add_carriage_return);
674 for (
auto const & [ref_name, ref_info] : views::
zip(header.ref_ids(), header.ref_id_info))
676 stream <<
"@SQ\tSN:";
680 stream <<
"\tLN:" << get<0>(ref_info);
682 if (!get<1>(ref_info).
empty())
683 stream <<
"\t" << get<1>(ref_info);
685 detail::write_eol(stream_it, options.add_carriage_return);
689 for (
auto const & read_group : header.read_groups)
692 <<
"\tID:" << get<0>(read_group);
694 if (!get<1>(read_group).
empty())
695 stream <<
"\t" << get<1>(read_group);
697 detail::write_eol(stream_it, options.add_carriage_return);
701 for (
auto const & program : header.program_infos)
704 <<
"\tID:" << program.id;
706 if (!program.name.empty())
707 stream <<
"\tPN:" << program.name;
709 if (!program.command_line_call.empty())
710 stream <<
"\tCL:" << program.command_line_call;
712 if (!program.previous.empty())
713 stream <<
"\tPP:" << program.previous;
715 if (!program.description.empty())
716 stream <<
"\tDS:" << program.description;
718 if (!program.version.empty())
719 stream <<
"\tVN:" << program.version;
721 detail::write_eol(stream_it, options.add_carriage_return);
725 for (
auto const &
comment : header.comments)
728 detail::write_eol(stream_it, options.add_carriage_return);
T back_inserter(T... args)
Provides seqan3::views::char_to.
Provides various utility functions.
constexpr auto assign_char_to
Assign a character to an alphabet object.
Definition alphabet/concept.hpp:521
@ comment
Comment field of arbitrary content, usually a string.
@ ref_id
The identifier of the (reference) sequence that seqan3::field::seq was aligned to.
constexpr auto is_char
Checks whether a given letter is the same as the template non-type argument.
Definition predicate.hpp:60
constexpr size_t size
The size of a type pack.
Definition type_pack/traits.hpp:143
seqan::stl::views::zip zip
A view adaptor that takes several views and returns tuple-like values from every i-th element of each...
Definition zip.hpp:24
Provides various utility functions.
Auxiliary functions for the SAM IO.
Provides seqan3::debug_stream and related types.
Provides seqan3::views::repeat_n.
Provides seqan3::views::slice.
Provides traits to inspect some information of a type, for example its name.
Provides seqan3::views::zip.