Provides files and formats for handling read mapping data. More...

Collaboration diagram for SAM File:

Classes
class	seqan3::format_bam
	The BAM format. More...

class	seqan3::format_sam
	The SAM format (tag). More...

struct	seqan3::ref_info_not_given
	Type tag which indicates that no reference information has been passed to the SAM file on construction. More...

class	seqan3::sam_file_header< ref_ids_type >
	Stores the header information of SAM/BAM files. More...

class	seqan3::sam_file_input< traits_type_, selected_field_ids_, valid_formats_ >
	A class for reading SAM files, both SAM and its binary representation BAM are supported. More...

struct	seqan3::sam_file_input_default_traits< ref_sequences_t, ref_ids_t >
	The default traits for seqan3::sam_file_input. More...

interface	seqan3::sam_file_input_format< t >
	The generic concept for alignment file input formats. More...

struct	seqan3::sam_file_input_options< sequence_legal_alphabet >
	The options type defines various option members that influence the behaviour of all or some formats. More...

interface	sam_file_input_traits
	The requirements a traits_type for seqan3::sam_file_input must meet. More...

class	seqan3::sam_file_output< selected_field_ids_, valid_formats_, ref_ids_type >
	A class for writing SAM files, both SAM and its binary representation BAM are supported. More...

interface	seqan3::sam_file_output_format< t >
	The generic concept for alignment file out formats. More...

struct	seqan3::sam_file_output_options
	The options type defines various option members that influence the behavior of all or some formats. More...

struct	seqan3::sam_file_program_info_t
	Stores information of the program/tool that was used to create a SAM/BAM file. More...

struct	seqan3::sam_flag_printer< sam_flag >
	A sam_flag can be printed as an integer value. More...

class	seqan3::sam_record< field_types, field_ids >
	The record type of seqan3::sam_file_input. More...

class	seqan3::sam_tag_dictionary
	The SAM tag dictionary class that stores all optional SAM fields. More...

struct	seqan3::sam_tag_type< tag_value >
	The generic base class. More...

Enumerations
enum class	seqan3::sam_flag : uint16_t { seqan3::sam_flag::none = 0 , seqan3::sam_flag::paired = 0x1 , seqan3::sam_flag::proper_pair = 0x2 , seqan3::sam_flag::unmapped = 0x4 , seqan3::sam_flag::mate_unmapped = 0x8 , seqan3::sam_flag::on_reverse_strand = 0x10 , seqan3::sam_flag::mate_on_reverse_strand = 0x20 , seqan3::sam_flag::first_in_pair = 0x40 , seqan3::sam_flag::second_in_pair = 0x80 , seqan3::sam_flag::secondary_alignment = 0x100 , seqan3::sam_flag::failed_filter = 0x200 , seqan3::sam_flag::duplicate = 0x400 , seqan3::sam_flag::supplementary_alignment = 0x800 }
	An enum flag that describes the properties of an aligned read (given as a SAM record). More...

Other literals
template<small_string< 2 > str>
constexpr uint16_t	seqan3::literals::operator""_tag ()
	The SAM tag literal, such that tags can be used in constant expressions.

template<small_string< 2 > str>
constexpr uint16_t	operator""_tag ()
	The SAM tag literal, such that tags can be used in constant expressions.

Detailed Description

Provides files and formats for handling read mapping data.

Introduction

SAM/BAM files are primarily used to store pairwise alignments of read mapping data.

Note: For a step-by-step guide take a look at our tutorial: SAM Input and Output in SeqAn.

The SAM file abstraction supports reading 10 different fields:

There exists one more field for SAM files, the seqan3::field::header_ptr, but this field is mostly used internally. Please see the seqan3::sam_file_output::header member function for details on how to access the seqan3::sam_file_header of the file.

All of these fields are retrieved by default (and in that order).

Reading SAM files

Construction and specialisation

The seqan3::sam_file_input class comes with four constructors: One for construction from a file name, one for construction from an existing stream and a known format and both of the former with or without additional reference information. Constructing from a file name automatically picks the format based on the extension of the file name. Constructing from a stream can be used if you have a non-file stream, like std::cin or std::istringstream. It also comes in handy, if you cannot use file-extension based detection, but know that your input file has a certain format.

Passing reference information, e.g.

ref_ids: The name of the references, e.g. "chr1", "chr2", ...
ref_sequences: The reference sequence information in the same order as the ref_ids.

comes in handy once you want to convert the CIGAR string, read from your file, into an actual alignment. This will be covered in the section "Transforming the CIGAR information into an actual alignment".In most cases the template parameters are deduced automatically:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <filesystem>
#include <sstream>
 
#include <seqan3/io/sam_file/input.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    // Create the temporary file.
    auto tmp_file = std::filesystem::temp_directory_path() / "my.sam";
    std::ofstream tmp_stream{tmp_file};
    tmp_stream << sam_file_raw;
    tmp_stream.close();
 
    seqan3::sam_file_input fin{tmp_file}; // SAM format assumed, regular std::ifstream taken as stream
    std::filesystem::remove(tmp_file);
}

Reading from a std::istringstream:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/io/sam_file/input.hpp>
 
auto input = R"(@HD VN:1.6  SO:coordinate
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{input}, seqan3::format_sam{}};
    //                          ^ no need to specify the template arguments
}

Note that this is not the same as writing sam_file_input<> (with angle brackets). In the latter case they are explicitly set to their default values, in the former case automatic deduction happens which chooses different parameters depending on the constructor arguments. For opening from file, sam_file_input<> would have also worked, but for opening from stream it would not have.

You can define your own traits type to further customise the types used by and returned by this class, see seqan3::sam_file_input_default_traits for more details. As mentioned above, specifying at least one template parameter yourself means that you loose automatic deduction. The following is equivalent to the automatic type deduction example with a stream from above:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto input = R"(@HD VN:1.6  SO:coordinate
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *)";
 
int main()
{
    // The default types; you can adjust this list if you don't want to read all this data.
    using default_fields = seqan3::fields<seqan3::field::seq,
                                          seqan3::field::id,
                                          seqan3::field::ref_id,
                                          seqan3::field::ref_offset,
                                          seqan3::field::cigar,
                                          seqan3::field::mapq,
                                          seqan3::field::qual,
                                          seqan3::field::flag,
                                          seqan3::field::mate,
                                          seqan3::field::tags,
                                          seqan3::field::header_ptr>;
 
    // The expected format:
    using sam_file_input_t = seqan3::sam_file_input<seqan3::sam_file_input_default_traits<>,
                                                    default_fields,
                                                    // Which formats are allowed:
                                                    seqan3::type_list<seqan3::format_sam>>;
 
    sam_file_input_t fin{std::istringstream{input}, seqan3::format_sam{}};
}

Reading record-wise

You can iterate over this file record-wise:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
 
    for (auto & rec : fin)
    {
        seqan3::debug_stream << "id:  " << rec.id() << '\n';
        seqan3::debug_stream << "read sequence: " << rec.sequence() << '\n';
        seqan3::debug_stream << "mapping position: " << rec.reference_position() << '\n';
        seqan3::debug_stream << "mapping quality: " << rec.mapping_quality() << '\n';
 
        // there are more fields read on default
    }
}

In the above example, rec has the type seqan3::sam_file_input::record_type which is a specialisation of seqan3::record and behaves like a std::tuple (that's why we can access it via get). Instead of using the seqan3::field based interface on the record, you could also use std::get<0> or even std::get<dna4_vector> to retrieve the sequence, but it is not recommended, because it is more error-prone.

Note: It is important to write auto & and not just auto, otherwise you will copy the record on every iteration. Since the buffer gets "refilled" on every iteration, you can also move the data out of the record if you want to store it somewhere without copying:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
#include <utility>
#include <vector>
 
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
 
    using record_type = typename decltype(fin)::record_type;
    std::vector<record_type> records{}; // store all my records in a vector
 
    for (auto & rec : fin)
        records.push_back(std::move(rec));
}

Reading record-wise (custom fields)

If you want to skip specific fields from the record you can pass a non-empty fields trait object to the seqan3::sam_file_input constructor to select the fields that should be read from the input. For example, you may only be interested in the mapping flag and mapping quality of your SAM data to get some statistics. The following snippets demonstrate the usage of such a fields trait object.

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw},
                               seqan3::format_sam{},
                               seqan3::fields<seqan3::field::flag, seqan3::field::mapq>{}};
 
    for (auto & rec : fin)
    {
        seqan3::debug_stream << "flag:  " << rec.flag() << '\n';
        seqan3::debug_stream << "mapping quality:  " << rec.mapping_quality() << '\n';
    }
}

When reading a file, all fields not present in the file (but requested implicitly or via the selected_field_ids parameter) are ignored and the respective value in the record stays empty.

Reading record-wise (decomposed records)

Instead of using get on the record, you can also use structured bindings to decompose the record into its elements. Considering the example of reading only the flag and mapping quality like before you can also write:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw},
                               seqan3::format_sam{},
                               seqan3::fields<seqan3::field::flag, seqan3::field::mapq>{}};
 
    for (auto & [flag, mapq] : fin) // the order is the same as specified in fields!
    {
        seqan3::debug_stream << "flag:  " << flag << '\n';
        seqan3::debug_stream << "mapping quality:  " << mapq << '\n';
    }
}

In this case you immediately get the two elements of the tuple: flag of seqan3::sam_file_input::flag_type and mapq of seqan3::sam_file_input::mapq_type.

Note: But beware: with structured bindings you do need to get the order of elements correctly!

Transforming the CIGAR information into an actual alignment

In SeqAn, we represent an alignment as a tuple of two seqan3::aligned_sequences.The conversion from a CIGAR string to an alignment can be done with the function seqan3::alignment_from_cigar. You need to pass the reference sequence with the position the read was aligned to and the read sequence. All of it is already in the record when reading a SAM file:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <seqan3/alignment/cigar_conversion/alignment_from_cigar.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/all.hpp>
 
using namespace seqan3::literals;
 
auto sam_file_raw = R"(@HD  VN:1.6
@SQ SN:ref  LN:34
read1   41  ref 1   61  1S1M1D1M1I  ref 10  300 ACGT    !##$    AS:i:2  NM:i:7
read2   42  ref 2   62  1H7M1D1M1S2H    ref 10  300 AGGCTGNAG   !##$&'()*   xy:B:S,3,4,5
read3   43  ref 3   63  1S1M1P1M1I1M1I1D1M1S    ref 10  300 GGAGTATA    !!*+,-./
)";
 
int main()
{
    // The reference sequence might be read from a different file.
    seqan3::dna5_vector reference = "ACTGATCGAGAGGATCTAGAGGAGATCGTAGGAC"_dna5;
 
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
    // You will probably read it from a file, e.g., like this:
    // seqan3::sam_file_input fin{"test.sam"};
 
    for (auto && rec : fin)
    {
        auto alignment =
            alignment_from_cigar(rec.cigar_sequence(), reference, rec.reference_position().value(), rec.sequence());
        seqan3::debug_stream << alignment << '\n';
    }
 
    // prints:
    // (ACT-,C-GT)
    // (CTGATCGAG,AGGCTGN-A)
    // (T-G-A-TC,G-AGTA-T)
}

The code will print the following:

(ACT-,C-GT)
(CTGATCGAG,AGGCTGN-A)
(T-G-A-TC,G-AGTA-T)

Views on files

Since SeqAn files are ranges, you can also create views over files. A useful example is to filter the records based on certain criteria, e.g. minimum length of the sequence field:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <ranges>
#include <sstream>
 
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/input.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
 
    auto minimum_length10_filter = std::views::filter(
        [](auto const & rec)
        {
            return std::ranges::size(rec.sequence()) >= 10;
        });
 
    for (auto & rec : fin | minimum_length10_filter) // only records with sequence length >= 10 will "appear"
        seqan3::debug_stream << rec.id() << '\n';
}

End of file

You can check whether a file is at its end by comparing begin() and end() (if they are the same, the file is at its end).

Formats

We currently support reading the following formats:

Writing SAM files

Construction and specialisation

The seqan3::sam_file_output class comes with two constructors, one for construction from a file name and one for construction from an existing stream and a known format. The first one automatically picks the format based on the extension of the file name. The second can be used if you have a non-file stream, like std::cout or std::ostringstream, that you want to read from and/or if you cannot use file-extension based detection, but know that your output file has a certain format.

In most cases the template parameters are deduced completely automatically:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <filesystem>
 
#include <seqan3/io/sam_file/output.hpp>
 
int main()
{
    auto tmp_file = std::filesystem::temp_directory_path() / "my.sam";
 
    seqan3::sam_file_output fout{tmp_file}; // SAM format detected, std::ofstream opened for file
 
    std::filesystem::remove(tmp_file);
}

Writing to std::cout:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <iostream>
 
#include <seqan3/io/sam_file/output.hpp>
 
int main()
{
    seqan3::sam_file_output fout{std::cout, seqan3::format_sam{}};
}

Note that this is not the same as writing sam_file_output<> (with angle brackets). In the latter case they are explicitly set to their default values, in the former case automatic deduction happens which chooses different parameters depending on the constructor arguments. For opening from file, sam_file_output<> would have also worked, but for opening from stream it would not have.

Writing record-wise

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
#include <string>
#include <vector>
 
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/io/sam_file/output.hpp>
#include <seqan3/utility/type_list/type_list.hpp>
 
int main()
{
    seqan3::sam_file_output fout{std::ostringstream{}, seqan3::format_sam{}};
 
    std::string ref_id;
    std::string read_id;
 
    std::vector<seqan3::dna5> read;
 
    // ... e.g. compute and alignment
 
    using alignment_type =
        std::pair<std::vector<seqan3::gapped<seqan3::dna5>>, std::vector<seqan3::gapped<seqan3::dna5>>>;
 
    alignment_type dummy_alignment{}; // an empty dummy alignment
 
    using types = seqan3::type_list<std::vector<seqan3::dna5>, std::string, alignment_type>;
    using types_as_ids = seqan3::fields<seqan3::field::seq, seqan3::field::id, seqan3::field::alignment>;
    // the record type specifies the fields we want to write
    using record_type = seqan3::record<types, types_as_ids>;
 
    // initialize record
    record_type rec{read, ref_id, dummy_alignment};
 
    // Write the record
    fout.push_back(rec);
 
    // same as
    fout.push_back(record_type{read, ref_id, dummy_alignment});
 
    // as all our fields are empty so this would print an
}

The easiest way to write to a SAM/BAM file is to use the push_back() member functions. These work similarly to how they work on a std::vector. You may also use a tuple like interface or the emplace_back() function but this is not recommended since one would have to keep track of the correct order of many fields (14 in total). For the record based interface using push_back() please also see the seqan3::record documentation on how to specify a record with the correct field and type lists.

You may also use the output file's iterator for writing, however, this rarely provides an advantage.

Writing record-wise (custom fields)

If you want to omit non-required parameter or change the order of the parameters, you can pass a non-empty fields trait object to the seqan3::sam_file_output constructor to select the fields that are used for interpreting the arguments.

The following snippet demonstrates the usage of such a field_traits object.

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <filesystem>
#include <sstream>
#include <tuple>
 
#include <seqan3/io/sam_file/output.hpp>
 
int main()
{
    // I only want to print the mapping position (field::ref_offset) and flag:
    seqan3::sam_file_output fout{std::ostringstream{},
                                 seqan3::format_sam{},
                                 seqan3::fields<seqan3::field::ref_offset, seqan3::field::flag>{}};
 
    unsigned mapping_pos{1300};
    seqan3::sam_flag flag{seqan3::sam_flag::none};
 
    // ...
 
    fout.emplace_back(mapping_pos, flag); // note that the order the arguments is now different, because
    // or:                                    you specified that REF_OFFSET should be first
    fout.push_back(std::tie(mapping_pos, flag));
}

A different way of passing custom fields to the file is to pass a seqan3::record – instead of a tuple – to push_back(). The seqan3::record clearly indicates which of its elements has which seqan3::field so the file will use that information instead of the template argument. This is especially handy when reading from one file and writing to another, because you don't have to configure the output file to match the input file, it will just work:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <filesystem>
#include <sstream>
 
#include <seqan3/io/sam_file/all.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  6H5M    *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    // fin uses custom fields, fout uses the default fields.
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw},
                               seqan3::format_sam{},
                               seqan3::fields<seqan3::field::ref_offset, seqan3::field::flag>{}};
    // output doesn't have to match the configuration of the input
    seqan3::sam_file_output fout{std::ostringstream{}, seqan3::format_sam{}};
 
    for (auto & r : fin)
        fout.push_back(r); // copy all the records.
}

This will copy the seqan3::field::flag and seqan3::field::ref_offset value into the new output file.

Note: Note that the other SAM columns in the output file will have a default value, so unless you specify to read all SAM columns (see seqan3::format_sam) the output file will not be equal to the input file.

Writing record-wise in batches

You can write multiple records at once, by assigning to the file:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
#include <string>
#include <tuple>
#include <vector>
 
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/io/sam_file/output.hpp>
 
int main()
{
    using namespace seqan3::literals;
 
    seqan3::sam_file_output fout{std::ostringstream{}, seqan3::format_sam{}};
 
    std::vector<std::tuple<seqan3::dna5_vector, std::string>> range{{"ACGT"_dna5, "First"},
                                                                    {"NATA"_dna5, "2nd"},
                                                                    {"GATA"_dna5, "Third"}}; // a range of "records"
 
    fout = range; // will iterate over the records and write them
    // equivalent to:
    range | fout;
}

File I/O pipelines

Record-wise writing in batches also works for writing from input files directly to output files, because input files are also input ranges in SeqAn:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <sstream>
 
#include <seqan3/io/sam_file/all.hpp>
 
auto sam_file_raw = R"(First    0   *   0   0   *   *   0   0   ACGT    *
2nd 0   *   0   0   *   *   0   0   NATA    *
Third   0   *   0   0   *   *   0   0   GATA    *
)";
 
int main()
{
    // copying a file in one line:
    seqan3::sam_file_output{std::ostringstream{}, seqan3::format_sam{}} =
        seqan3::sam_file_input{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
 
    // with seqan3::sam_file_output as a variable:
    seqan3::sam_file_output fout{std::ostringstream{}, seqan3::format_sam{}};
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
    fout = fin;
 
    // or in pipe notation:
    seqan3::sam_file_input{std::istringstream{sam_file_raw}, seqan3::format_sam{}}
        | seqan3::sam_file_output{std::ostringstream{}, seqan3::format_sam{}};
}

This can be combined with file-based views to create I/O pipelines:

#include <ranges>
#include <sstream>
 
#include <seqan3/io/sam_file/all.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  *   =   37  39  TTAGATAAAGGATACTG   *
r003    0   ref 29  30  *   *   0   0   GCCTAAGCTAA *   SA:Z:ref,29,-,6H5M,17,0;
r003    2064    ref 29  17  *   *   0   0   TAGGC   *   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  *   =   7   -39 CAGCGGCAT   *   NM:i:1
)";
 
int main()
{
    auto input_file = seqan3::sam_file_input{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
    input_file | std::views::take(3) // take only the first 3 records
        | seqan3::sam_file_output{std::ostringstream{}, seqan3::format_sam{}};
}

Formats

We currently support writing the following formats:

Enumeration Type Documentation

◆ sam_flag

enum class seqan3::sam_flag : uint16_t

strong

An enum flag that describes the properties of an aligned read (given as a SAM record).

See also: seqan3::enum_bitwise_operators enables combining enum values.

The SAM flag are bitwise flags, which means that each value corresponds to a specific bit that is set and that they can be combined and tested using binary operations. See this tutorial for an introduction on bitwise operations on enum flags.

Example:

// SPDX-FileCopyrightText: 2006-2024 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2024 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: CC0-1.0
 
#include <iostream>
#include <sstream>
 
#include <seqan3/io/sam_file/all.hpp>
 
auto sam_file_raw = R"(@HD  VN:1.6  SO:coordinate   GO:none
@SQ SN:ref  LN:45
r001    99  ref 7   30  8M2I4M1D3M  =   37  39  TTAGATAAAGGATACTG   !!!!!!!!!!!!!!!!!
r003    0   ref 29  30  5S6M    *   0   0   GCCTAAGCTAA !!!!!!!!!!! SA:Z:ref,29,-,6H5M,17,0;
r003    4   *   29  17  *   *   0   0   TAGGC   @@@@@   SA:Z:ref,9,+,5S6M,30,1;
r001    147 ref 237 30  9M  =   7   -39 CAGCGGCAT   !!!!!!!!!   NM:i:1
)";
 
int main()
{
    seqan3::sam_file_input fin{std::istringstream{sam_file_raw}, seqan3::format_sam{}};
 
    for (auto & rec : fin)
    {
        // Check if a certain flag value (bit) is set:
        if (static_cast<bool>(rec.flag() & seqan3::sam_flag::unmapped))
            std::cout << "Read " << rec.id() << " is unmapped\n";
 
        if (rec.base_qualities()[0] < seqan3::assign_char_to('@', seqan3::phred42{})) // low quality
        {
            // Set a flag value (bit):
            rec.flag() |= seqan3::sam_flag::failed_filter;
            // Note that this does not affect other flag values (bits),
            // e.g. `rec.flag() & seqan3::sam_flag::unmapped` may still be true
        }
 
        // Unset a flag value (bit):
        rec.flag() &= ~seqan3::sam_flag::duplicate; // not marked as a duplicate anymore
    }
}

Adapted from the SAM specifications are the following additional information to some flag values:

For each read/contig in a SAM file, it is required that one and only one line associated with the read has neither the seqan3::sam_flag::secondary_alignment nor the seqan3::sam_flag::supplementary_alignment flag value set (satisfies FLAG & 0x900 == 0). This line is called the primary alignment of the read.
seqan3::sam_flag::secondary_alignment (bit 0x100) marks the alignment not to be used in certain analyses when the tools in use are aware of this bit. It is typically used to flag alternative mappings when multiple mappings are presented in a SAM.
seqan3::sam_flag::supplementary_alignment (bit 0x800) indicates that the corresponding alignment line is part of a chimeric alignment. If the SAM/BAM file corresponds to long reads (nanopore/pacbio) this happens when reads are split before being aligned and the best matching part is marked as primary, while all other aligned parts are marked supplementary.
seqan3::sam_flag::unmapped (bit 0x4) is the only reliable place to tell whether the read is unmapped. If seqan3::sam_flag::unmapped is set, no assumptions can be made about RNAME, POS, CIGAR, MAPQ, and seqan3::sam_flag::proper_pair, seqan3::sam_flag::secondary_alignment, and seqan3::sam_flag::supplementary_alignment (bits 0x2, 0x100, and 0x800).
seqan3::sam_flag::on_reverse_strand (bit 0x10) indicates whether the read sequence has been reverse complemented and the quality string is reversed. When bit seqan3::sam_flag::unmapped (0x4) is unset, this corresponds to the strand to which the segment has been mapped: seqan3::sam_flag::on_reverse_strand (bit 0x10) unset indicates the forward strand, while set indicates the reverse strand. When seqan3::sam_flag::unmapped (0x4) is set, this indicates whether the unmapped read is stored in its original orientation as it came off the sequencing machine.
seqan3::sam_flag::first_in_pair and seqan3::sam_flag::second_in_pair (bits 0x40 and 0x80) reflect the read ordering within each template inherent in the sequencing technology used. If seqan3::sam_flag::first_in_pair and seqan3::sam_flag::second_in_pair (0x40 and 0x80) are both set, the read is part of a linear template, but it is neither the first nor the last read. If both are unset, the index of the read in the template is unknown. This may happen for a non-linear template or when this information is lost during data processing.
If seqan3::sam_flag::paired (bit 0x1) is unset, no assumptions can be made about seqan3::sam_flag::proper_pair, seqan3::sam_flag::mate_unmapped, seqan3::sam_flag::mate_on_reverse_strand, seqan3::sam_flag::first_in_pair and seqan3::sam_flag::second_in_pair (bits 0x2, 0x8, 0x20, 0x40 and 0x80).

See also: https://broadinstitute.github.io/picard/explain-flags.html

Remarks: For a complete overview, take a look at SAM File

Enumerator
none	None of the flags below are set.
paired	The aligned read is paired (paired-end sequencing).
proper_pair	The two aligned reads in a pair have a proper distance between each other.
unmapped	The read is not mapped to a reference (unaligned).
mate_unmapped	The mate of this read is not mapped to a reference (unaligned).
on_reverse_strand	The read sequence has been reverse complemented before being mapped (aligned).
mate_on_reverse_strand	The mate sequence has been reverse complemented before being mapped (aligned).
first_in_pair	Indicates the ordering (see details in the seqan3::sam_flag description).
second_in_pair	Indicates the ordering (see details in the seqan3::sam_flag description).
secondary_alignment	This read alignment is an alternative (possibly suboptimal) to the primary.
failed_filter	The read alignment failed a filter, e.g. quality controls.
duplicate	The read is marked as a PCR duplicate or optical duplicate.
supplementary_alignment	This sequence is part of a split alignment and is not the primary alignment.

Function Documentation

◆ operator""_tag() [1/2]

template<small_string< 2 > str>

constexpr uint16_t seqan3::literals::operator""_tag ( )

constexpr

The SAM tag literal, such that tags can be used in constant expressions.

Template Parameters

char_t The char type. Usually char. Parameter pack ...s must be of length 2 since SAM tags consist of two letters (char0 and char1).

Returns: The unique identifier of the SAM tag computed by char0 * 128 + char1.

A SAM tag consists of two letters, initialized via the string literal ""_tag, which delegate to its unique id.

#include <seqan3/io/sam_file/sam_tag_dictionary.hpp>
 
using namespace seqan3::literals;
 
// ...
 
uint16_t tag_id = "NM"_tag; // tag_id = 10061

The purpose of those tags is to fill or query the seqan3::sam_tag_dictionary for a specific key (tag_id) and retrieve the corresponding value.

◆ operator""_tag() [2/2]

template<small_string< 2 > str>

constexpr uint16_t operator""_tag ( )

Template Parameters

char_t The char type. Usually char. Parameter pack ...s must be of length 2 since SAM tags consist of two letters (char0 and char1).

Returns: The unique identifier of the SAM tag computed by char0 * 128 + char1.

A SAM tag consists of two letters, initialized via the string literal ""_tag, which delegate to its unique id.

#include <seqan3/io/sam_file/sam_tag_dictionary.hpp>
 
using namespace seqan3::literals;
 
// ...
 
uint16_t tag_id = "NM"_tag; // tag_id = 10061

The purpose of those tags is to fill or query the seqan3::sam_tag_dictionary for a specific key (tag_id) and retrieve the corresponding value.

Classes

Enumerations

Other literals

Detailed Description

Introduction

Reading SAM files

Construction and specialisation

Reading record-wise

Reading record-wise (custom fields)

Reading record-wise (decomposed records)

Transforming the CIGAR information into an actual alignment

Views on files

End of file

Formats

Writing SAM files

Construction and specialisation

Writing record-wise

Writing record-wise (custom fields)

Writing record-wise in batches

File I/O pipelines

Formats

Enumeration Type Documentation

◆ sam_flag

Function Documentation

◆ operator""_tag() [1/2]

◆ operator""_tag() [2/2]