SeqAn3 3.4.0-rc.4
The Modern C++ library for sequence analysis.
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
format_sam.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2006-2025 Knut Reinert & Freie Universität Berlin
2// SPDX-FileCopyrightText: 2016-2025 Knut Reinert & MPI für molekulare Genetik
3// SPDX-License-Identifier: BSD-3-Clause
4
10#pragma once
11
12#include <iterator>
13#include <ranges>
14#include <string>
15#include <vector>
16
17#include <seqan3/alphabet/concept.hpp>
18#include <seqan3/alphabet/views/to_char.hpp>
19#include <seqan3/io/sam_file/detail/cigar.hpp>
20#include <seqan3/io/sam_file/detail/format_sam_base.hpp>
21#include <seqan3/io/sam_file/header.hpp>
22#include <seqan3/io/sam_file/input_format_concept.hpp>
23#include <seqan3/io/sam_file/output_options.hpp>
24#include <seqan3/io/sam_file/sam_flag.hpp>
25#include <seqan3/io/sam_file/sam_tag_dictionary.hpp>
26#include <seqan3/io/sequence_file/input_format_concept.hpp>
27#include <seqan3/io/sequence_file/output_options.hpp>
28#include <seqan3/io/stream/detail/fast_ostreambuf_iterator.hpp>
29#include <seqan3/io/views/detail/istreambuf_view.hpp>
30#include <seqan3/io/views/detail/take_until_view.hpp>
31#include <seqan3/utility/detail/type_name_as_string.hpp>
32#include <seqan3/utility/range/to.hpp>
33#include <seqan3/utility/tuple/concept.hpp>
34#include <seqan3/utility/views/slice.hpp>
35
36namespace seqan3
37{
38
104class format_sam : protected detail::format_sam_base
105{
106public:
110 // construction cannot be noexcept because this class has a std::string variable as a quality string buffer.
111 format_sam() = default;
112 format_sam(format_sam const &) = delete;
113 format_sam & operator=(format_sam const &) = delete;
114 format_sam(format_sam &&) = default;
115 format_sam & operator=(format_sam &&) = default;
116 ~format_sam() = default;
117
119
121 static inline std::vector<std::string> file_extensions{
122 {"sam"},
123 };
124
125protected:
126 template <typename stream_type, // constraints checked by file
127 typename seq_legal_alph_type,
128 typename stream_pos_type,
129 typename seq_type, // other constraints checked inside function
130 typename id_type,
131 typename qual_type>
132 void read_sequence_record(stream_type & stream,
133 sequence_file_input_options<seq_legal_alph_type> const & options,
134 stream_pos_type & position_buffer,
135 seq_type & sequence,
136 id_type & id,
137 qual_type & qualities);
138
139 template <typename stream_type, // constraints checked by file
140 typename seq_type, // other constraints checked inside function
141 typename id_type,
142 typename qual_type>
143 void write_sequence_record(stream_type & stream,
144 sequence_file_output_options const & SEQAN3_DOXYGEN_ONLY(options),
145 seq_type && sequence,
146 id_type && id,
147 qual_type && qualities);
148
149 template <typename stream_type, // constraints checked by file
150 typename seq_legal_alph_type,
151 typename ref_seqs_type,
152 typename ref_ids_type,
153 typename stream_pos_type,
154 typename seq_type,
155 typename id_type,
156 typename ref_seq_type,
157 typename ref_id_type,
158 typename ref_offset_type,
159 typename cigar_type,
160 typename flag_type,
161 typename mapq_type,
162 typename qual_type,
163 typename mate_type,
164 typename tag_dict_type,
165 typename e_value_type,
166 typename bit_score_type>
167 void read_alignment_record(stream_type & stream,
168 sam_file_input_options<seq_legal_alph_type> const & SEQAN3_DOXYGEN_ONLY(options),
169 ref_seqs_type & ref_seqs,
170 sam_file_header<ref_ids_type> & header,
171 stream_pos_type & position_buffer,
172 seq_type & seq,
173 qual_type & qual,
174 id_type & id,
175 ref_seq_type & SEQAN3_DOXYGEN_ONLY(ref_seq),
176 ref_id_type & ref_id,
177 ref_offset_type & ref_offset,
178 cigar_type & cigar_vector,
179 flag_type & flag,
180 mapq_type & mapq,
181 mate_type & mate,
182 tag_dict_type & tag_dict,
183 e_value_type & SEQAN3_DOXYGEN_ONLY(e_value),
184 bit_score_type & SEQAN3_DOXYGEN_ONLY(bit_score));
185
186 template <typename stream_type,
187 typename header_type,
188 typename seq_type,
189 typename id_type,
190 typename ref_seq_type,
191 typename ref_id_type,
192 typename qual_type,
193 typename mate_type,
194 typename tag_dict_type,
195 typename e_value_type,
196 typename bit_score_type>
197 void write_alignment_record(stream_type & stream,
198 sam_file_output_options const & options,
199 header_type && header,
200 seq_type && seq,
201 qual_type && qual,
202 id_type && id,
203 ref_seq_type && SEQAN3_DOXYGEN_ONLY(ref_seq),
204 ref_id_type && ref_id,
205 std::optional<int32_t> ref_offset,
206 std::vector<cigar> const & cigar_vector,
207 sam_flag const flag,
208 uint8_t const mapq,
209 mate_type && mate,
210 tag_dict_type && tag_dict,
211 e_value_type && SEQAN3_DOXYGEN_ONLY(e_value),
212 bit_score_type && SEQAN3_DOXYGEN_ONLY(bit_score));
213
214private:
216 std::string tmp_qual{};
217
219 static constexpr std::string_view dummy{};
220
222 sam_file_header<> default_header{};
223
225 std::array<std::string_view, 11> raw_record{};
226
228 std::string_view const & default_or(detail::ignore_t) const noexcept
229 {
230 return dummy;
231 }
232
234 template <typename t>
235 decltype(auto) default_or(t && v) const noexcept
236 {
237 return std::forward<t>(v);
238 }
239
240 template <arithmetic value_type>
241 void read_sam_dict_vector(seqan3::detail::sam_tag_variant & variant, std::string_view const str, value_type value);
242
243 void read_sam_byte_vector(seqan3::detail::sam_tag_variant & variant, std::string_view const str);
244
245 void read_sam_dict(std::string_view const tag_str, sam_tag_dictionary & target);
246
247 template <typename stream_it_t, std::ranges::forward_range field_type>
248 void write_range_or_asterisk(stream_it_t & stream_it, field_type && field_value);
249
250 template <typename stream_it_t>
251 void write_range_or_asterisk(stream_it_t & stream_it, char const * const field_value);
252
253 template <typename stream_it_t>
254 void write_tag_fields(stream_it_t & stream, sam_tag_dictionary const & tag_dict, char const separator);
255};
256
258template <typename stream_type, // constraints checked by file
259 typename seq_legal_alph_type,
260 typename stream_pos_type,
261 typename seq_type, // other constraints checked inside function
262 typename id_type,
263 typename qual_type>
264inline void format_sam::read_sequence_record(stream_type & stream,
265 sequence_file_input_options<seq_legal_alph_type> const & options,
266 stream_pos_type & position_buffer,
267 seq_type & sequence,
268 id_type & id,
269 qual_type & qualities)
270{
271 sam_file_input_options<seq_legal_alph_type> align_options;
272
273 {
274 read_alignment_record(stream,
275 align_options,
276 std::ignore,
277 default_header,
278 position_buffer,
279 sequence,
280 qualities,
281 id,
282 std::ignore,
283 std::ignore,
284 std::ignore,
285 std::ignore,
286 std::ignore,
287 std::ignore,
288 std::ignore,
289 std::ignore,
290 std::ignore,
291 std::ignore);
292 }
293
294 if constexpr (!detail::decays_to_ignore_v<seq_type>)
295 if (std::ranges::distance(sequence) == 0)
296 throw parse_error{"The sequence information must not be empty."};
297 if constexpr (!detail::decays_to_ignore_v<id_type>)
298 {
299 if (std::ranges::distance(id) == 0)
300 throw parse_error{"The id information must not be empty."};
301 if (options.truncate_ids)
302 id = id | detail::take_until_and_consume(is_space) | ranges::to<id_type>();
303 }
304}
305
307template <typename stream_type, // constraints checked by file
308 typename seq_type, // other constraints checked inside function
309 typename id_type,
310 typename qual_type>
311inline void format_sam::write_sequence_record(stream_type & stream,
312 sequence_file_output_options const & SEQAN3_DOXYGEN_ONLY(options),
313 seq_type && sequence,
314 id_type && id,
315 qual_type && qualities)
316{
317 using default_mate_t = std::tuple<std::string_view, std::optional<int32_t>, int32_t>;
318
319 sam_file_output_options output_options;
320
321 write_alignment_record(stream,
322 output_options,
323 /*header*/ std::ignore,
324 /*seq*/ default_or(sequence),
325 /*qual*/ default_or(qualities),
326 /*id*/ default_or(id),
327 /*ref_seq*/ std::string_view{},
328 /*ref_id*/ std::string_view{},
329 /*ref_offset*/ -1,
330 /*cigar_vector*/ std::vector<cigar>{},
331 /*flag*/ sam_flag::none,
332 /*mapq*/ 0,
333 /*mate*/ default_mate_t{},
334 /*tag_dict*/ sam_tag_dictionary{},
335 /*e_value*/ 0,
336 /*bit_score*/ 0);
337}
338
340template <typename stream_type, // constraints checked by file
341 typename seq_legal_alph_type,
342 typename ref_seqs_type,
343 typename ref_ids_type,
344 typename stream_pos_type,
345 typename seq_type,
346 typename id_type,
347 typename ref_seq_type,
348 typename ref_id_type,
349 typename ref_offset_type,
350 typename cigar_type,
351 typename flag_type,
352 typename mapq_type,
353 typename qual_type,
354 typename mate_type,
355 typename tag_dict_type,
356 typename e_value_type,
357 typename bit_score_type>
358inline void
359format_sam::read_alignment_record(stream_type & stream,
360 sam_file_input_options<seq_legal_alph_type> const & SEQAN3_DOXYGEN_ONLY(options),
361 ref_seqs_type & ref_seqs,
362 sam_file_header<ref_ids_type> & header,
363 stream_pos_type & position_buffer,
364 seq_type & seq,
365 qual_type & qual,
366 id_type & id,
367 ref_seq_type & SEQAN3_DOXYGEN_ONLY(ref_seq),
368 ref_id_type & ref_id,
369 ref_offset_type & ref_offset,
370 cigar_type & cigar_vector,
371 flag_type & flag,
372 mapq_type & mapq,
373 mate_type & mate,
374 tag_dict_type & tag_dict,
375 e_value_type & SEQAN3_DOXYGEN_ONLY(e_value),
376 bit_score_type & SEQAN3_DOXYGEN_ONLY(bit_score))
377{
378 static_assert(detail::decays_to_ignore_v<ref_offset_type>
379 || detail::is_type_specialisation_of_v<ref_offset_type, std::optional>,
380 "The ref_offset must be a specialisation of std::optional.");
381
382 auto stream_it = detail::fast_istreambuf_iterator{*stream.rdbuf()};
383
384 auto stream_view = detail::istreambuf(stream);
385
386 int32_t ref_offset_tmp{}; // needed to read the ref_offset (int) beofre storing it in std::optional<uint32_t>
387 std::ranges::range_value_t<decltype(header.ref_ids())> ref_id_tmp{}; // in case mate is requested but ref_offset not
388
389 // Header
390 // -------------------------------------------------------------------------------------------------------------
391 if (is_char<'@'>(*stream_it)) // we always read the header if present
392 {
393 read_header(stream_view, header, ref_seqs);
394
395 if (std::ranges::begin(stream_view) == std::ranges::end(stream_view)) // file has no records
396 return;
397 }
398
399 // Store the current file position in the buffer.
400 position_buffer = stream.tellg();
401
402 // We don't know wether we have 11 or 12 fields, since the tags are optional.
403 // The last field will thus contain either the quality sequence
404 // or the quality sequence AND tags. This will be handled at the respective fields below.
405 stream_it.cache_record_into('\n', '\t', raw_record);
406
407 // Fields 1-5: ID FLAG REF_ID REF_OFFSET MAPQ
408 // -------------------------------------------------------------------------------------------------------------
409 if constexpr (!detail::decays_to_ignore_v<id_type>)
410 read_forward_range_field(raw_record[0], id);
411
412 uint16_t flag_integral{};
413 read_arithmetic_field(raw_record[1], flag_integral);
414 flag = sam_flag{flag_integral};
415
416 read_forward_range_field(raw_record[2], ref_id_tmp);
417 check_and_assign_ref_id(ref_id, ref_id_tmp, header, ref_seqs);
418
419 read_arithmetic_field(raw_record[3], ref_offset_tmp);
420 --ref_offset_tmp; // SAM format is 1-based but SeqAn operates 0-based
421
422 if (ref_offset_tmp == -1)
423 ref_offset = std::nullopt; // indicates an unmapped read -> ref_offset is not set
424 else if (ref_offset_tmp > -1)
425 ref_offset = ref_offset_tmp;
426 else if (ref_offset_tmp < -1)
427 throw format_error{"No negative values are allowed for field::ref_offset."};
428
429 if constexpr (!detail::decays_to_ignore_v<mapq_type>)
430 read_arithmetic_field(raw_record[4], mapq);
431
432 // Field 6: CIGAR
433 // -------------------------------------------------------------------------------------------------------------
434 if constexpr (!detail::decays_to_ignore_v<cigar_type>)
435 cigar_vector = detail::parse_cigar(raw_record[5]);
436
437 // Field 7-9: (RNEXT PNEXT TLEN) = MATE
438 // -------------------------------------------------------------------------------------------------------------
439 if constexpr (!detail::decays_to_ignore_v<mate_type>)
440 {
441 std::ranges::range_value_t<decltype(header.ref_ids())> tmp_mate_ref_id{};
442 read_forward_range_field(raw_record[6], tmp_mate_ref_id); // RNEXT
443
444 if (tmp_mate_ref_id == "=") // indicates "same as ref id"
445 {
446 if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
447 get<0>(mate) = ref_id;
448 else
449 check_and_assign_ref_id(get<0>(mate), ref_id_tmp, header, ref_seqs);
450 }
451 else
452 {
453 check_and_assign_ref_id(get<0>(mate), tmp_mate_ref_id, header, ref_seqs);
454 }
455
456 int32_t tmp_pnext{};
457 read_arithmetic_field(raw_record[7], tmp_pnext); // PNEXT
458
459 if (tmp_pnext > 0)
460 get<1>(mate) = --tmp_pnext; // SAM format is 1-based but SeqAn operates 0-based.
461 else if (tmp_pnext < 0)
462 throw format_error{"No negative values are allowed at the mate mapping position."};
463 // tmp_pnext == 0 indicates an unmapped mate -> do not fill std::optional get<1>(mate)
464
465 read_arithmetic_field(raw_record[8], get<2>(mate)); // TLEN
466 }
467
468 // Field 10: Sequence
469 // -------------------------------------------------------------------------------------------------------------
470 if constexpr (!detail::decays_to_ignore_v<seq_type>)
471 {
472 std::string_view const seq_str = raw_record[9];
473
474 if (!seq_str.starts_with('*')) // * indicates missing sequence information
475 {
476 seq.resize(seq_str.size());
477 constexpr auto is_legal_alph = char_is_valid_for<seq_legal_alph_type>;
478
479 for (size_t i = 0; i < seq_str.size(); ++i)
480 {
481 if (!is_legal_alph(seq_str[i]))
482 throw parse_error{std::string{"Encountered an unexpected letter: "} + "char_is_valid_for<"
483 + detail::type_name_as_string<seq_legal_alph_type>
484 + "> evaluated to false on " + detail::make_printable(seq_str[i])};
485
486 seq[i] = assign_char_to(seq_str[i], std::ranges::range_value_t<seq_type>{});
487 }
488 }
489 }
490
491 // Field 11: Quality
492 // -------------------------------------------------------------------------------------------------------------
493 // We don't know wether we have 11 or 12 fields, since the tags are optional.
494 // The last field will thus contain either the quality sequence
495 // or the quality sequence AND tags.
496 size_t tag_begin_pos = raw_record[10].find('\t');
497
498 std::string_view qualities =
499 (tag_begin_pos == std::string_view::npos) ? raw_record[10] : raw_record[10].substr(0, tag_begin_pos);
500
501 if constexpr (!detail::decays_to_ignore_v<qual_type>)
502 read_forward_range_field(qualities, qual);
503
504 if constexpr (!detail::decays_to_ignore_v<seq_type> && !detail::decays_to_ignore_v<qual_type>)
505 {
506 if (std::ranges::distance(seq) != 0 && std::ranges::distance(qual) != 0
507 && std::ranges::distance(seq) != std::ranges::distance(qual))
508 {
509 throw format_error{detail::to_string("Sequence length (",
510 std::ranges::distance(seq),
511 ") and quality length (",
512 std::ranges::distance(qual),
513 ") must be the same.")};
514 }
515 }
516
517 // All remaining optional fields if any: SAM tags dictionary
518 // -------------------------------------------------------------------------------------------------------------
519 if constexpr (!detail::decays_to_ignore_v<tag_dict_type>)
520 {
521 while (tag_begin_pos != std::string_view::npos) // read all tags if present
522 {
523 ++tag_begin_pos; // skip '\t'
524 size_t const tag_end_pos = raw_record[10].find('\t', tag_begin_pos);
525
526 char const * tag_begin = raw_record[10].begin() + tag_begin_pos;
527 char const * tag_end =
528 (tag_end_pos == std::string_view::npos) ? raw_record[10].end() : raw_record[10].begin() + tag_end_pos;
529
530 read_sam_dict(std::string_view{tag_begin, tag_end}, tag_dict);
531
532 tag_begin_pos = tag_end_pos;
533 }
534 }
535
536 assert(stream_it == std::default_sentinel_t{} || *stream_it == '\n');
537 ++stream_it; // Move from end of record to the beginning of the next or to the end of the stream.
538}
539
541template <typename stream_type,
542 typename header_type,
543 typename seq_type,
544 typename id_type,
545 typename ref_seq_type,
546 typename ref_id_type,
547 typename qual_type,
548 typename mate_type,
549 typename tag_dict_type,
550 typename e_value_type,
551 typename bit_score_type>
552inline void format_sam::write_alignment_record(stream_type & stream,
553 sam_file_output_options const & options,
554 header_type && header,
555 seq_type && seq,
556 qual_type && qual,
557 id_type && id,
558 ref_seq_type && SEQAN3_DOXYGEN_ONLY(ref_seq),
559 ref_id_type && ref_id,
560 std::optional<int32_t> ref_offset,
561 std::vector<cigar> const & cigar_vector,
562 sam_flag const flag,
563 uint8_t const mapq,
564 mate_type && mate,
565 tag_dict_type && tag_dict,
566 e_value_type && SEQAN3_DOXYGEN_ONLY(e_value),
567 bit_score_type && SEQAN3_DOXYGEN_ONLY(bit_score))
568{
569 /* Note the following general things:
570 *
571 * - Given the SAM specifications, all fields may be empty
572 *
573 * - arithmetic values default to 0 while all others default to '*'
574 *
575 * - Because of the former, arithmetic values can be directly streamed
576 * into 'stream' as operator<< is defined for all arithmetic types
577 * and the default value (0) is also the SAM default.
578 *
579 * - All other non-arithmetic values need to be checked for emptiness
580 */
581
582 // ---------------------------------------------------------------------
583 // Type Requirements (as static asserts for user friendliness)
584 // ---------------------------------------------------------------------
585 static_assert((std::ranges::forward_range<seq_type> && alphabet<std::ranges::range_reference_t<seq_type>>),
586 "The seq object must be a std::ranges::forward_range over "
587 "letters that model seqan3::alphabet.");
588
589 static_assert((std::ranges::forward_range<id_type> && alphabet<std::ranges::range_reference_t<id_type>>),
590 "The id object must be a std::ranges::forward_range over "
591 "letters that model seqan3::alphabet.");
592
593 if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
594 {
595 static_assert((std::ranges::forward_range<ref_id_type> || std::integral<std::remove_reference_t<ref_id_type>>
596 || detail::is_type_specialisation_of_v<std::remove_cvref_t<ref_id_type>, std::optional>),
597 "The ref_id object must be a std::ranges::forward_range "
598 "over letters that model seqan3::alphabet.");
599
600 if constexpr (std::integral<std::remove_cvref_t<ref_id_type>>
601 || detail::is_type_specialisation_of_v<std::remove_cvref_t<ref_id_type>, std::optional>)
602 static_assert(!detail::decays_to_ignore_v<header_type>,
603 "If you give indices as reference id information the header must also be present.");
604 }
605
606 static_assert((std::ranges::forward_range<qual_type> && alphabet<std::ranges::range_reference_t<qual_type>>),
607 "The qual object must be a std::ranges::forward_range "
608 "over letters that model seqan3::alphabet.");
609
610 static_assert(tuple_like<std::remove_cvref_t<mate_type>>,
611 "The mate object must be a std::tuple of size 3 with "
612 "1) a std::ranges::forward_range with a value_type modelling seqan3::alphabet, "
613 "2) a std::integral or std::optional<std::integral>, and "
614 "3) a std::integral.");
615
616 static_assert(
617 ((std::ranges::forward_range<decltype(std::get<0>(mate))>
618 || std::integral<std::remove_cvref_t<decltype(std::get<0>(mate))>>
619 || detail::is_type_specialisation_of_v<std::remove_cvref_t<decltype(std::get<0>(mate))>, std::optional>)
620 && (std::integral<std::remove_cvref_t<decltype(std::get<1>(mate))>>
621 || detail::is_type_specialisation_of_v<std::remove_cvref_t<decltype(std::get<1>(mate))>, std::optional>)
622 && std::integral<std::remove_cvref_t<decltype(std::get<2>(mate))>>),
623 "The mate object must be a std::tuple of size 3 with "
624 "1) a std::ranges::forward_range with a value_type modelling seqan3::alphabet, "
625 "2) a std::integral or std::optional<std::integral>, and "
626 "3) a std::integral.");
627
628 if constexpr (std::integral<std::remove_cvref_t<decltype(std::get<0>(mate))>>
629 || detail::is_type_specialisation_of_v<std::remove_cvref_t<decltype(std::get<0>(mate))>,
630 std::optional>)
631 static_assert(!detail::decays_to_ignore_v<header_type>,
632 "If you give indices as mate reference id information the header must also be present.");
633
634 static_assert(std::same_as<std::remove_cvref_t<tag_dict_type>, sam_tag_dictionary>,
635 "The tag_dict object must be of type seqan3::sam_tag_dictionary.");
636
637 // ---------------------------------------------------------------------
638 // logical Requirements
639 // ---------------------------------------------------------------------
640 if constexpr (!detail::decays_to_ignore_v<header_type> && !detail::decays_to_ignore_v<ref_id_type>
641 && !std::integral<std::remove_reference_t<ref_id_type>>
642 && !detail::is_type_specialisation_of_v<std::remove_reference_t<ref_id_type>, std::optional>)
643 {
644
645 if (options.sam_require_header && !std::ranges::empty(ref_id))
646 {
647 auto id_it = header.ref_dict.end();
648
649 if constexpr (std::ranges::contiguous_range<decltype(ref_id)> && std::ranges::sized_range<decltype(ref_id)>
650 && std::ranges::borrowed_range<decltype(ref_id)>)
651 {
652 id_it = header.ref_dict.find(std::span{std::ranges::data(ref_id), std::ranges::size(ref_id)});
653 }
654 else
655 {
656 using header_ref_id_type = std::remove_reference_t<decltype(header.ref_ids()[0])>;
657
658 static_assert(implicitly_convertible_to<ref_id_type, header_ref_id_type>,
659 "The ref_id type is not convertible to the reference id information stored in the "
660 "reference dictionary of the header object.");
661
662 id_it = header.ref_dict.find(ref_id);
663 }
664
665 if (id_it == header.ref_dict.end()) // no reference id matched
666 throw format_error{detail::to_string("The ref_id '",
667 ref_id,
668 "' was not in the list of references:",
669 header.ref_ids())};
670 }
671 }
672
673 if (ref_offset.has_value() && (ref_offset.value() + 1) < 0)
674 throw format_error{"The ref_offset object must be a std::integral >= 0."};
675
676 // ---------------------------------------------------------------------
677 // Writing the Header on first call
678 // ---------------------------------------------------------------------
679 if constexpr (!detail::decays_to_ignore_v<header_type>)
680 {
681 if (options.sam_require_header && !header_was_written)
682 {
683 write_header(stream, options, header);
684 header_was_written = true;
685 }
686 }
687
688 // ---------------------------------------------------------------------
689 // Writing the Record
690 // ---------------------------------------------------------------------
691
692 detail::fast_ostreambuf_iterator stream_it{*stream.rdbuf()};
693 constexpr char separator{'\t'};
694
695 write_range_or_asterisk(stream_it, id);
696 *stream_it = separator;
697
698 stream_it.write_number(static_cast<uint16_t>(flag));
699 *stream_it = separator;
700
701 if constexpr (!detail::decays_to_ignore_v<ref_id_type>)
702 {
703 if constexpr (std::integral<std::remove_reference_t<ref_id_type>>)
704 {
705 write_range_or_asterisk(stream_it, (header.ref_ids())[ref_id]);
706 }
707 else if constexpr (detail::is_type_specialisation_of_v<std::remove_reference_t<ref_id_type>, std::optional>)
708 {
709 if (ref_id.has_value())
710 write_range_or_asterisk(stream_it, (header.ref_ids())[ref_id.value()]);
711 else
712 *stream_it = '*';
713 }
714 else
715 {
716 write_range_or_asterisk(stream_it, ref_id);
717 }
718 }
719 else
720 {
721 *stream_it = '*';
722 }
723
724 *stream_it = separator;
725
726 // SAM is 1 based, 0 indicates unmapped read if optional is not set
727 stream_it.write_number(ref_offset.value_or(-1) + 1);
728 *stream_it = separator;
729
730 stream_it.write_number(static_cast<unsigned>(mapq));
731 *stream_it = separator;
732
733 if (!std::ranges::empty(cigar_vector))
734 {
735 for (auto & c : cigar_vector) //TODO THIS IS PROBABLY TERRIBLE PERFORMANCE_WISE
736 stream_it.write_range(c.to_string());
737 }
738 else
739 {
740 *stream_it = '*';
741 }
742
743 *stream_it = separator;
744
745 if constexpr (std::integral<std::remove_reference_t<decltype(get<0>(mate))>>)
746 {
747 write_range_or_asterisk(stream_it, (header.ref_ids())[get<0>(mate)]);
748 }
749 else if constexpr (detail::is_type_specialisation_of_v<std::remove_reference_t<decltype(get<0>(mate))>,
750 std::optional>)
751 {
752 if (get<0>(mate).has_value())
753 write_range_or_asterisk(stream_it, header.ref_ids()[get<0>(mate).value()]);
754 else
755 *stream_it = '*';
756 }
757 else
758 {
759 write_range_or_asterisk(stream_it, get<0>(mate));
760 }
761
762 *stream_it = separator;
763
764 if constexpr (detail::is_type_specialisation_of_v<std::remove_cvref_t<decltype(get<1>(mate))>, std::optional>)
765 {
766 // SAM is 1 based, 0 indicates unmapped read if optional is not set
767 stream_it.write_number(get<1>(mate).value_or(-1) + 1);
768 *stream_it = separator;
769 }
770 else
771 {
772 stream_it.write_number(get<1>(mate));
773 *stream_it = separator;
774 }
775
776 stream_it.write_number(get<2>(mate));
777 *stream_it = separator;
778
779 write_range_or_asterisk(stream_it, seq);
780 *stream_it = separator;
781
782 write_range_or_asterisk(stream_it, qual);
783
784 write_tag_fields(stream_it, tag_dict, separator);
785
786 stream_it.write_end_of_line(options.add_carriage_return);
787}
788
805template <arithmetic value_type>
806inline void format_sam::read_sam_dict_vector(seqan3::detail::sam_tag_variant & variant,
807 std::string_view const str,
808 value_type value)
809{
810 std::vector<value_type> tmp_vector{};
811 size_t start_pos{0};
812 size_t end_pos{0};
813
814 while (start_pos != std::string_view::npos)
815 {
816 end_pos = str.find(',', start_pos);
817 auto end = (end_pos == std::string_view::npos) ? str.end() : str.begin() + end_pos;
818 read_arithmetic_field(std::string_view{str.begin() + start_pos, end}, value);
819 tmp_vector.push_back(value);
820
821 start_pos = (end_pos == std::string_view::npos) ? end_pos : end_pos + 1;
822 }
823 variant = std::move(tmp_vector);
824}
825
837inline void format_sam::read_sam_byte_vector(seqan3::detail::sam_tag_variant & variant, std::string_view const str)
838{
839 std::vector<std::byte> tmp_vector{};
840 // std::from_chars cannot directly parse into a std::byte
841 uint8_t dummy_byte{};
842
843 if (str.size() % 2 != 0)
844 throw format_error{"[CORRUPTED SAM FILE] Hexadecimal tag must have even number of digits."};
845
846 // H encodes bytes in a hexadecimal format. Two hex values are stored for each byte as characters.
847 // E.g., '1' and 'A' need one byte each and are read as `\x1A`, which is 27 in decimal.
848 for (auto hex_begin = str.begin(), hex_end = str.begin() + 2; hex_begin != str.end(); hex_begin += 2, hex_end += 2)
849 {
850 auto res = std::from_chars(hex_begin, hex_end, dummy_byte, 16);
851
852 if (res.ec == std::errc::invalid_argument)
853 throw format_error{std::string("[CORRUPTED SAM FILE] The string '") + std::string(hex_begin, hex_end)
854 + "' could not be cast into type uint8_t."};
855
856 if (res.ec == std::errc::result_out_of_range)
857 throw format_error{std::string("[CORRUPTED SAM FILE] Casting '") + std::string(str)
858 + "' into type uint8_t would cause an overflow."};
859
860 tmp_vector.push_back(std::byte{dummy_byte});
861 }
862
863 variant = std::move(tmp_vector);
864}
865
881inline void format_sam::read_sam_dict(std::string_view const tag_str, sam_tag_dictionary & target)
882{
883 /* Every SAM tag has the format "[TAG]:[TYPE_ID]:[VALUE]", where TAG is a two letter
884 name tag which is converted to a unique integer identifier and TYPE_ID is one character in [A,i,Z,H,B,f]
885 describing the type for the upcoming VALUES. If TYPE_ID=='B' it signals an array of comma separated
886 VALUE's and the inner value type is identified by the character following ':', one of [cCsSiIf].
887 */
888 assert(tag_str.size() > 5);
889
890 uint16_t tag = static_cast<uint16_t>(tag_str[0]) << 8;
891 tag += static_cast<uint16_t>(tag_str[1]);
892
893 char type_id = tag_str[3];
894
895 switch (type_id)
896 {
897 case 'A': // char
898 {
899 assert(tag_str.size() == 6);
900 target[tag] = tag_str[5];
901 break;
902 }
903 case 'i': // int32_t
904 {
905 int32_t tmp;
906 read_arithmetic_field(tag_str.substr(5), tmp);
907 target[tag] = tmp;
908 break;
909 }
910 case 'f': // float
911 {
912 float tmp;
913 read_arithmetic_field(tag_str.substr(5), tmp);
914 target[tag] = tmp;
915 break;
916 }
917 case 'Z': // string
918 {
919 target[tag] = std::string{tag_str.substr(5)};
920 break;
921 }
922 case 'H':
923 {
924 read_sam_byte_vector(target[tag], tag_str.substr(5));
925 break;
926 }
927 case 'B': // Array. Value type depends on second char [cCsSiIf]
928 {
929 assert(tag_str.size() > 6);
930 char array_value_type_id = tag_str[5];
931
932 switch (array_value_type_id)
933 {
934 case 'c': // int8_t
935 read_sam_dict_vector(target[tag], tag_str.substr(7), int8_t{});
936 break;
937 case 'C': // uint8_t
938 read_sam_dict_vector(target[tag], tag_str.substr(7), uint8_t{});
939 break;
940 case 's': // int16_t
941 read_sam_dict_vector(target[tag], tag_str.substr(7), int16_t{});
942 break;
943 case 'S': // uint16_t
944 read_sam_dict_vector(target[tag], tag_str.substr(7), uint16_t{});
945 break;
946 case 'i': // int32_t
947 read_sam_dict_vector(target[tag], tag_str.substr(7), int32_t{});
948 break;
949 case 'I': // uint32_t
950 read_sam_dict_vector(target[tag], tag_str.substr(7), uint32_t{});
951 break;
952 case 'f': // float
953 read_sam_dict_vector(target[tag], tag_str.substr(7), float{});
954 break;
955 default:
956 throw format_error{std::string("The first character in the numerical ")
957 + "id of a SAM tag must be one of [cCsSiIf] but '" + array_value_type_id
958 + "' was given."};
959 }
960 break;
961 }
962 default:
963 throw format_error{std::string("The second character in the numerical id of a "
964 "SAM tag ([TAG]:[TYPE_ID]:[VALUE]) must be one of [A,i,Z,H,B,f] but '")
965 + type_id + "' was given."};
966 }
967}
968
976template <typename stream_it_t, std::ranges::forward_range field_type>
977inline void format_sam::write_range_or_asterisk(stream_it_t & stream_it, field_type && field_value)
978{
979 if (std::ranges::empty(field_value))
980 {
981 *stream_it = '*';
982 }
983 else
984 {
985 if constexpr (std::same_as<std::remove_cvref_t<std::ranges::range_reference_t<field_type>>, char>)
986 stream_it.write_range(field_value);
987 else // convert from alphabets to their character representation
988 stream_it.write_range(field_value | views::to_char);
989 }
990}
991
998template <typename stream_it_t>
999inline void format_sam::write_range_or_asterisk(stream_it_t & stream_it, char const * const field_value)
1000{
1001 write_range_or_asterisk(stream_it, std::string_view{field_value});
1002}
1003
1011template <typename stream_it_t>
1012inline void
1013format_sam::write_tag_fields(stream_it_t & stream_it, sam_tag_dictionary const & tag_dict, char const separator)
1014{
1015 auto const stream_variant_fn = [&stream_it](auto && arg) // helper to print a std::variant
1016 {
1017 using T = std::remove_cvref_t<decltype(arg)>;
1018
1019 if constexpr (std::ranges::input_range<T>)
1020 {
1021 if constexpr (std::same_as<std::remove_cvref_t<std::ranges::range_reference_t<T>>, char>)
1022 {
1023 stream_it.write_range(arg);
1024 }
1025 else if constexpr (std::same_as<std::remove_cvref_t<std::ranges::range_reference_t<T>>, std::byte>)
1026 {
1027 if (!std::ranges::empty(arg))
1028 {
1029 stream_it.write_number(std::to_integer<uint8_t>(*std::ranges::begin(arg)));
1030
1031 for (auto && elem : arg | std::views::drop(1))
1032 {
1033 *stream_it = ',';
1034 stream_it.write_number(std::to_integer<uint8_t>(elem));
1035 }
1036 }
1037 }
1038 else
1039 {
1040 if (!std::ranges::empty(arg))
1041 {
1042 stream_it.write_number(*std::ranges::begin(arg));
1043
1044 for (auto && elem : arg | std::views::drop(1))
1045 {
1046 *stream_it = ',';
1047 stream_it.write_number(elem);
1048 }
1049 }
1050 }
1051 }
1052 else if constexpr (std::same_as<std::remove_cvref_t<T>, char>)
1053 {
1054 *stream_it = arg;
1055 }
1056 else // number
1057 {
1058 stream_it.write_number(arg);
1059 }
1060 };
1061
1062 for (auto & [tag, variant] : tag_dict)
1063 {
1064 *stream_it = separator;
1065
1066 char const char0 = tag / 256;
1067 char const char1 = tag % 256;
1068
1069 *stream_it = char0;
1070 *stream_it = char1;
1071 *stream_it = ':';
1072 *stream_it = detail::sam_tag_type_char[variant.index()];
1073 *stream_it = ':';
1074
1075 if (detail::sam_tag_type_char_extra[variant.index()] != '\0')
1076 {
1077 *stream_it = detail::sam_tag_type_char_extra[variant.index()];
1078 *stream_it = ',';
1079 }
1080
1081 std::visit(stream_variant_fn, variant);
1082 }
1083}
1084
1085} // namespace seqan3
concept.hpp
Core alphabet concept and free function/type trait wrappers.
std::ranges::begin
T begin(T... args)
seqan3::debug_stream_type
A "pretty printer" for most SeqAn data structures and related types.
Definition debug_stream_type.hpp:79
seqan3::format_sam
The SAM format (tag).
Definition format_sam.hpp:105
seqan3::format_sam::~format_sam
~format_sam()=default
Defaulted.
seqan3::format_sam::operator=
format_sam & operator=(format_sam const &)=delete
Deleted. Header holds a unique_ptr.
seqan3::format_sam::file_extensions
static std::vector< std::string > file_extensions
The valid file extensions for this format; note that you can modify this value.
Definition format_sam.hpp:121
seqan3::format_sam::write_sequence_record
void write_sequence_record(stream_type &stream, sequence_file_output_options const &options, seq_type &&sequence, id_type &&id, qual_type &&qualities)
Write the given fields to the specified stream.
Definition format_sam.hpp:311
seqan3::format_sam::write_alignment_record
void write_alignment_record(stream_type &stream, sam_file_output_options const &options, header_type &&header, seq_type &&seq, qual_type &&qual, id_type &&id, ref_seq_type &&ref_seq, ref_id_type &&ref_id, std::optional< int32_t > ref_offset, std::vector< cigar > const &cigar_vector, sam_flag const flag, uint8_t const mapq, mate_type &&mate, tag_dict_type &&tag_dict, e_value_type &&e_value, bit_score_type &&bit_score)
Write the given fields to the specified stream.
Definition format_sam.hpp:552
seqan3::format_sam::format_sam
format_sam(format_sam &&)=default
Defaulted.
seqan3::format_sam::operator=
format_sam & operator=(format_sam &&)=default
Defaulted.
seqan3::format_sam::read_alignment_record
void read_alignment_record(stream_type &stream, sam_file_input_options< seq_legal_alph_type > const &options, ref_seqs_type &ref_seqs, sam_file_header< ref_ids_type > &header, stream_pos_type &position_buffer, seq_type &seq, qual_type &qual, id_type &id, ref_seq_type &ref_seq, ref_id_type &ref_id, ref_offset_type &ref_offset, cigar_type &cigar_vector, flag_type &flag, mapq_type &mapq, mate_type &mate, tag_dict_type &tag_dict, e_value_type &e_value, bit_score_type &bit_score)
Read from the specified stream and back-insert into the given field buffers.
Definition format_sam.hpp:359
seqan3::format_sam::read_sequence_record
void read_sequence_record(stream_type &stream, sequence_file_input_options< seq_legal_alph_type > const &options, stream_pos_type &position_buffer, seq_type &sequence, id_type &id, qual_type &qualities)
Read from the specified stream and back-insert into the given field buffers.
Definition format_sam.hpp:264
seqan3::format_sam::format_sam
format_sam(format_sam const &)=delete
Deleted. Header holds a unique_ptr.
seqan3::format_sam::format_sam
format_sam()=default
Defaulted.
seqan3::sam_tag_dictionary
The SAM tag dictionary class that stores all optional SAM fields.
Definition sam_tag_dictionary.hpp:327
std::end
T end(T... args)
fast_ostreambuf_iterator.hpp
Provides seqan3::detail::fast_ostreambuf_iterator.
std::string_view::find
T find(T... args)
format_sam_base.hpp
Provides the seqan3::format_sam_base that can be inherited from.
std::from_chars
T from_chars(T... args)
seqan3::views::to_char
auto const to_char
A view that calls seqan3::to_char() on each element in the input range.
Definition to_char.hpp:60
seqan3::assign_char_to
constexpr auto assign_char_to
Assign a character to an alphabet object.
Definition alphabet/concept.hpp:517
seqan3::sam_flag
sam_flag
An enum flag that describes the properties of an aligned read (given as a SAM record).
Definition sam_flag.hpp:73
seqan3::sam_flag::none
@ none
None of the flags below are set.
seqan3::field::flag
@ flag
The alignment flag (bit information), uint16_t value.
seqan3::field::ref_offset
@ ref_offset
Sequence (seqan3::field::ref_seq) relative start position (0-based), unsigned value.
seqan3::field::ref_seq
@ ref_seq
The (reference) "sequence" information, usually a range of nucleotides or amino acids.
seqan3::field::mapq
@ mapq
The mapping quality of the seqan3::field::seq alignment, usually a Phred-scaled score.
seqan3::field::bit_score
@ bit_score
The bit score (statistical significance indicator), unsigned value.
seqan3::field::mate
@ mate
The mate pair information given as a std::tuple of reference name, offset and template length.
seqan3::field::ref_id
@ ref_id
The identifier of the (reference) sequence that seqan3::field::seq was aligned to.
seqan3::field::seq
@ seq
The "sequence", usually a range of nucleotides or amino acids.
seqan3::field::qual
@ qual
The qualities, usually in Phred score notation.
seqan3::is_char
constexpr auto is_char
Checks whether a given letter is the same as the template non-type argument.
Definition predicate.hpp:60
seqan3::is_space
constexpr auto is_space
Checks whether c is a space character.
Definition predicate.hpp:122
seqan3::ranges::to
seqan::stl::ranges::to to
Converts a range to a container. <dl class="no-api">This entity is not part of the SeqAn API....
Definition to.hpp:23
seqan3::list_traits::drop
typename decltype(detail::split_after< i >(list_t{}))::second_type drop
Return a seqan3::type_list of the types in the input type list, except the first n.
Definition type_list/traits.hpp:392
header.hpp
Provides the seqan3::sam_file_header class.
std::variant::index
T index(T... args)
sequence
The generic concept for a (biological) sequence.
cigar.hpp
Auxiliary functions for the SAM IO.
istreambuf_view.hpp
Provides seqan3::detail::istreambuf.
seqan3
The main SeqAn3 namespace.
Definition aligned_sequence_concept.hpp:26
std
SeqAn specific customisations in the standard namespace.
std::string::push_back
T push_back(T... args)
input_format_concept.hpp
Provides seqan3::sam_file_input_format and auxiliary classes.
output_options.hpp
Provides seqan3::sam_file_output_options.
sam_flag.hpp
Provides helper data structures for the seqan3::sam_file_output.
sam_tag_dictionary.hpp
Provides the seqan3::sam_tag_dictionary class and auxiliaries.
input_format_concept.hpp
Provides seqan3::sequence_file_input_format and auxiliary classes.
output_options.hpp
Provides seqan3::sequence_file_output_options.
std::string_view::size
T size(T... args)
slice.hpp
Provides seqan3::views::slice.
seqan3::format_error
Thrown if information given to output format didn't match expectations.
Definition io/exception.hpp:88
seqan3::parse_error
Thrown if there is a parse error, such as reading an unexpected character from an input stream.
Definition io/exception.hpp:45
seqan3::sam_file_output_options
The options type defines various option members that influence the behavior of all or some formats.
Definition sam_file/output_options.hpp:23
seqan3::sam_file_output_options::add_carriage_return
bool add_carriage_return
The default plain text line-ending is "\n", but on Windows an additional carriage return is recommend...
Definition sam_file/output_options.hpp:27
seqan3::sam_file_output_options::sam_require_header
bool sam_require_header
Whether to require a header for SAM files.
Definition sam_file/output_options.hpp:41
seqan3::sequence_file_output_options
The options type defines various option members that influence the behaviour of all or some formats.
Definition sequence_file/output_options.hpp:23
std::string_view::substr
T substr(T... args)
take_until_view.hpp
Provides seqan3::views::take_until and seqan3::views::take_until_or_throw.
to.hpp
Provides seqan3::ranges::to.
to_char.hpp
Provides seqan3::views::to_char.
type_name_as_string.hpp
Provides traits to inspect some information of a type, for example its name.
concept.hpp
Provides seqan3::tuple_like.
std::visit
T visit(T... args)
Hide me