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