SeqAn3  3.0.2
The Modern C++ library for sequence analysis.
Alignment Input and Output in SeqAn

Learning Objective:

Learning Objective:
You will get an overview of how to read and write alignment files. This tutorial is a walk-through with links into the API documentation and also meant as a source for copy-and-paste code.

DifficultyMedium
Duration60 min
Prerequisite tutorialsQuick Setup (using CMake), Alphabets in SeqAn, Sequence File Input and Output
Recommended reading

Introduction

Alignment files are used to store pairwise alignments between two (biological) sequences. Common file formats are the Sequence Alignment/Map format (SAM) and BLAST output format. Next to the alignment, those formats store additional information like the start positions or mapping qualities. Alignment files are a little more complex than sequence files but the basic design is the same. If you are new to SeqAn, we strongly recommend to do the tutorial Sequence File Input and Output first.

Alignment file formats

SAM format

SAM stands for Sequence Alignment/Map format. It is a TAB-delimited text format consisting of a header section, which is optional, and an alignment section (see the official SAM specifications).

Here is an example of a SAM file:

//!\file
@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA *
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 5M * 0 0 TAGGC *
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1

The following table summarises the columns of a SAM file:

# SAM Column ID Description
1 QNAME Query template NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
4 POS 1-based leftmost mapping POSition
5 MAPQ MAPping Quality
6 CIGAR CIGAR string
7 RNEXT Reference name of the mate/next read
8 PNEXT Position of the mate/next read
9 TLEN observed Template LENgth
10 SEQ segment SEQuence
11 QUAL ASCII of Phred-scaled base QUALity+33

If you want to read more about the SAM format, take a look at the official specifications.

BAM format

BAM is the binary format version of SAM. It provides the same data as the SAM format with the only difference that the header is mandatory.

Alignment file fields

The Alignment file abstraction supports reading 15 different fields:

  1. seqan3::field::seq
  2. seqan3::field::id
  3. seqan3::field::offset
  4. seqan3::field::ref_seq
  5. seqan3::field::ref_id
  6. seqan3::field::ref_offset
  7. seqan3::field::alignment
  8. seqan3::field::cigar
  9. seqan3::field::mapq
  10. seqan3::field::qual
  11. seqan3::field::flag
  12. seqan3::field::mate
  13. seqan3::field::tags
  14. seqan3::field::evalue
  15. seqan3::field::bit_score

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

Note that some of the fields are specific to the SAM format, while some are specific to BLAST. To make things clearer, here is the table of SAM columns connected to the corresponding alignment file field:

# SAM Column ID FIELD name
1 QNAME seqan3::field::id
2 FLAG seqan3::field::flag
3 RNAME seqan3::field::ref_id
4 POS seqan3::field::ref_offset
5 MAPQ seqan3::field::mapq
6 CIGAR implicitly stored in seqan3::field::alignment or directly in seqan3::field::cigar
7 RNEXT seqan3::field::mate (tuple pos 0)
8 PNEXT seqan3::field::mate (tuple pos 1)
9 TLEN seqan3::field::mate (tuple pos 2)
10 SEQ seqan3::field::seq
11 QUAL seqan3::field::qual

File extensions

The formerly introduced formats can be identified by the following file name extensions (this is important for automatic format detection from a file name as you will learn in the next section).

File Format File Extensions
SAM .sam
BAM .bam

You can access and modify the valid file extensions via the file_extension member variable in a format tag:

Reading alignment files

Before we start, you should copy and paste this example file into a file location of your choice (we use /tmp/ in the examples, so make sure you adjust your path).

Attention
Make sure the file you copied is tab delimited!

Construction

The construction works analogously to sequence files by passing a file name, in which case all template parameters are automatically deduced (by the file name extension). Or you can pass a stream (e.g. std::cin or std::stringstream), but then you need to know your format beforehand:

int main()
{
seqan3::alignment_file_input fin_from_filename{"/tmp/my.sam"};
}

Reading custom fields

In many cases you are not interested in all of the information in a file. For this purpose, we provide the possibility to select specific seqan3::field's for a file. The file will read only those fields and fill the record accordingly.

You can select fields by providing a seqan3::fields object as an extra parameter to the constructor:

int main()
{
auto filename = std::filesystem::temp_directory_path()/"example.sam";
for (auto & [id, seq, flag /*order!*/] : fin)
{
seqan3::debug_stream << id << '\n';
seqan3::debug_stream << flag << '\n';
}
}
Attention
The order in which you specify the selected fields determines the order of elements in the seqan3::record.

In the example above we only select the id, sequence and flag information so the seqan3::record object has three tuple elements that are decomposed using structural bindings.

Note that this is possible for all SeqAn file objects.

Assignment 1: Accumulating mapping qualities

Let's assume we want to compute the average mapping quality of a SAM file.

For this purpose, write a small program that

  • only reads the mapping quality (field::mapq) out of a SAM file and
  • computes the average of all qualities.

Use the following file to test your program:

@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1

It should output:

Average: 27.4

Solution

int main()
{
std::filesystem::path tmp_dir = std::filesystem::temp_directory_path(); // get the temp directory
double sum{};
size_t c{};
std::ranges::for_each(fin.begin(), fin.end(), [&sum, &c] (auto & rec)
{
sum += seqan3::get<seqan3::field::mapq>(rec);
++c;
});
seqan3::debug_stream << "Average: " << (sum/c) << '\n';
}

Alignment representation in the SAM format

In SeqAn we represent an alignment as a tuple of two aligned_sequences, as you have probably learned by now from the alignment tutorial.

The SAM format is the common output format of read mappers where you align short read sequences to one or more large reference sequences. In fact, the SAM format stores those alignment information only partially: It does not store the reference sequence but only the read sequence and a CIGAR string representing the alignment based on the read.

Take this SAM record as an example:

r003 73 ref 3 17 1M1D4M * 0 0 TAGGC *

The record gives you the following information: A read with name r003 has been mapped to a reference with name ref at position 3 (in the reference, counting from 1) with a quality of 17 (Phred scaled). The flag has a value of 73 which indicates that the read is paired, the first in pair, but the mate is unmapped (see this website for a nice explanation of SAM flags). Fields set to 0 or * are defaulted and contain no information.

The cigar string is 1M1D4M which represents the following alignment:

1 2 3 4 5 6 7 8 9 ...
ref N N N N N N N N N ...
read T - A G G C

where the reference sequence is not known (represented by N). You will learn in the next section how to handle additional reference sequence information.

If you want to read up more about cigar strings, take a look at the SAM specifications or the SAMtools paper.

Reading the CIGAR string

By default, the seqan3::alignment_file_input will always read the seqan3::field::cigar and store it into a std::vector<seqan3::cigar>:

int main()
{
std::filesystem::path tmp_dir = std::filesystem::temp_directory_path(); // get the temp directory
seqan3::alignment_file_input fin{tmp_dir/"my.sam"}; // default fields
for (auto & rec : fin)
seqan3::debug_stream << seqan3::get<seqan3::field::cigar>(rec) << '\n'; // access cigar vector
}

Reading the CIGAR information into an actual alignment

In SeqAn, the conversion from a CIGAR string to an alignment (two aligned_sequences) is done automatically for you. You can access it by querying seqan3::field::alignment from the record:

int main()
{
auto filename = std::filesystem::temp_directory_path()/"example.sam";
for (auto & [ id, alignment ] : fin)
{
seqan3::debug_stream << id << ": " << std::get<1>(alignment) << '\n';
}
}

In the example above, you can only safely access the aligned read.

Attention
The unknown aligned reference sequence at the first position in the alignment tuple cannot be accessed (e.g. via the operator[]). It is represented by a dummy type that throws on access.

Although the SAM format does not handle reference sequence information, you can provide these information to the seqan3::alignment_file_input which automatically fills the alignment object. You can pass reference ids and reference sequences as additional constructor parameters:

int main()
{
using seqan3::operator""_dna5;
auto filename = std::filesystem::temp_directory_path()/"example.sam";
std::vector<std::string> ref_ids{"ref"}; // list of one reference name
std::vector<seqan3::dna5_vector> ref_sequences{"AGAGTTCGAGATCGAGGACTAGCGACGAGGCAGCGAGCGATCGAT"_dna5};
for (auto & [ alignment ] : fin)
{
seqan3::debug_stream << alignment << '\n'; // Now you can print the whole alignment!
}
}

Assignment 2: Combining sequence and alignment files

Read in the following reference sequence FASTA file (see the sequence file tutorial if you need a reminder):

>chr1
ACAGCAGGCATCTATCGGCGGATCGATCAGGCAGGCAGCTACTGG
>chr2
ACAGCAGGCATCTATCGGCGGATCGATCAGGCAGGCAGCTACTGTAATGGCATCAAAATCGGCATG

Then read in the following SAM file while providing the reference sequence information. Only read in the id, reference id, mapping quality, and alignment.

@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:45
@SQ SN:chr2 LN:66
r001 99 chr1 7 60 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r003 0 chr1 9 60 5S6M * 0 0 GCCTAAGCTAA *
r004 0 chr2 16 60 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 chr2 18 10 5M * 0 0 TAGGC *

With those information do the following:

  • Filter the alignment records and only take those with a mapping quality >= 30. (Take a look at the tutorial Reading paired-end reads for a reminder how to use views on files)
  • For the resulting alignments, print which read was mapped against with reference id and the number of seqan3::gap's in each sequence (aligned reference and read sequence).
Note
reference ids (field::ref_id) are given as an index of type std::optional<int32_t> that denote the position of the reference id in the ref_ids vector passed to the alignment file.

Your program should print the following:

r001 mapped against 0 with 1 gaps in the read sequence and 2 gaps in the reference sequence.
r003 mapped against 0 with 0 gaps in the read sequence and 0 gaps in the reference sequence.
r004 mapped against 1 with 14 gaps in the read sequence and 0 gaps in the reference sequence.

Solution

#include <string>
#include <vector>
int main()
{
std::filesystem::path tmp_dir = std::filesystem::temp_directory_path(); // get the temp directory
// read in reference information
seqan3::sequence_file_input reference_file{tmp_dir/"reference.fasta"};
for (auto && record : reference_file)
{
ref_ids.push_back(std::move(seqan3::get<seqan3::field::id>(record)));
ref_seqs.push_back(std::move(seqan3::get<seqan3::field::seq>(record)));
}
using field_type = seqan3::fields<seqan3::field::id,
seqan3::field::ref_id,
seqan3::alignment_file_input mapping_file{tmp_dir/"mapping.sam", ref_ids, ref_seqs, field_type{}};
#if !SEQAN3_WORKAROUND_GCC_93983
auto mapq_filter = std::views::filter([] (auto & rec) { return seqan3::get<seqan3::field::mapq>(rec) >= 30; });
#endif // !SEQAN3_WORKAROUND_GCC_93983
#if SEQAN3_WORKAROUND_GCC_93983
for (auto & [id, ref_id, mapq, alignment] : mapping_file /*| mapq_filter*/)
#else // ^^^ workaround / no workaround vvv
for (auto & [id, ref_id, mapq, alignment] : mapping_file | mapq_filter)
#endif // SEQAN3_WORKAROUND_GCC_93983
{
using seqan3::get;
size_t sum_ref{};
for (auto const & char_ref : get<0>(alignment))
if (char_ref == seqan3::gap{})
++sum_ref;
size_t sum_read{};
for (auto const & char_read : get<1>(alignment))
if (char_read == seqan3::gap{})
++sum_read;
seqan3::debug_stream << id << " mapped against " << ref_id << " with "
<< sum_read << " gaps in the read sequence and "
<< sum_ref << " gaps in the reference sequence.\n";
}
}

Writing alignment files

Writing records

When writing a SAM file without any further specifications, the default file assumes that all fields are provided. Since those are quite a lot for alignment files, we usually want to write only a subset of the data stored in the SAM format and default the rest.

For this purpose, you can also select specific fields by giving an additional seqan3::fields object to the constructor.

Attention
The order of the field tags in your seqan3::fields object will determine the order of values stored in the record type!
int main()
{
auto filename = std::filesystem::temp_directory_path()/"out.sam";
size_t mymapq{0};
seqan3::sam_flag flag{seqan3::sam_flag::unmapped};
// ...
fout.emplace_back(flag, mymapq);
// or:
fout.push_back(std::tie(flag, mymapq));
}

Note that this only works because in the SAM format all fields are optional. So if we provide less fields when writing, default values are printed.

Assignment 3: Writing id and sequence information

Write a small program that writes the following read ids + sequences:

read1: ACGATCGACTAGCTACGATCAGCTAGCAG
read2: AGAAAGAGCGAGGCTATTTTAGCGAGTTA

Your ids can be of type std::string and your sequences of type std::vector<seqan3::dna4>.

Your resulting SAM file should look like this:

read1 0 * 0 0 * * 0 0 ACGATCGACTAGCTACGATCAGCTAGCAG *
read2 0 * 0 0 * * 0 0 AGAAAGAGCGAGGCTATTTTAGCGAGTTA *

Solution

using seqan3::operator""_dna4;
int main()
{
std::vector<std::string> ids = {"read1", "read2"};
std::vector<std::vector<seqan3::dna4>> seqs = {"ACGATCGACTAGCTACGATCAGCTAGCAG"_dna4,
"AGAAAGAGCGAGGCTATTTTAGCGAGTTA"_dna4};
for (size_t i = 0; i < ids.size(); ++i)
{
fout.emplace_back(ids[i], seqs[i]);
}
}

debug_stream.hpp
Provides seqan3::debug_stream and related types.
seqan3::gap
The alphabet of a gap character '-'.
Definition: gap.hpp:37
dna4.hpp
Provides seqan3::dna4, container aliases and string literals.
std::filesystem::temp_directory_path
T temp_directory_path(T... args)
std::vector< std::string >
std::vector::size
T size(T... args)
seqan3::get
constexpr auto const & get(configuration< configs_t... > const &config) noexcept
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition: configuration.hpp:627
seqan3::field::id
@ id
The identifier, usually a string.
seqan3::alignment_file_output
A class for writing alignment files, e.g. SAM, BAL, BLAST, ...
Definition: output.hpp:163
seqan3::sam_flag
sam_flag
An enum flag that describes the properties of an aligned read (given as a SAM record).
Definition: misc.hpp:71
record.hpp
Provides the seqan3::record template and the seqan3::field enum.
input.hpp
Provides seqan3::sequence_file_input and corresponding traits classes.
filesystem
This header includes C++17 filesystem support and imports it into namespace std::filesystem (independ...
std::filesystem::path
seqan3::fields
A class template that holds a choice of seqan3::field.
Definition: record.hpp:166
std::tie
T tie(T... args)
std::vector::push_back
T push_back(T... args)
seqan3::seq
constexpr sequenced_policy seq
Global execution policy object for sequenced execution policy.
Definition: execution.hpp:54
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:142
seqan3::alignment_file_input
A class for reading alignment files, e.g. SAM, BAM, BLAST ...
Definition: input.hpp:353
seqan3::views::move
auto const move
A view that turns lvalue-references into rvalue-references.
Definition: move.hpp:68
seqan3::field::mapq
@ mapq
The mapping quality of the SEQ alignment, usually a ohred-scaled score.
dna5.hpp
Provides seqan3::dna5, container aliases and string literals.
gap.hpp
Provides seqan3::gap.
alignment_file_input_traits::ref_sequences
using ref_sequences
The type of range over reference sequences; must model std::ranges::forward_range,...
alignment_file_input_traits::ref_ids
using ref_ids
The type of range over reference sequences; must model std::ranges::forward_range,...
ranges
Adaptations of concepts from the Ranges TS.
seqan3::sequence_file_input
A class for reading sequence files, e.g. FASTA, FASTQ ...
Definition: input.hpp:316
input.hpp
Provides seqan3::alignment_file_input and corresponding traits classes.
seqan3::format_sam
The SAM format (tag).
Definition: format_sam.hpp:126
seqan3::debug_stream
debug_stream_type debug_stream
A global instance of seqan3::debug_stream_type.
Definition: debug_stream.hpp:42
seqan3::field::flag
@ flag
The alignment flag (bit information), uint16_t value.
all.hpp
Meta-include for the alignment IO submodule.
std::cin
std::ref
T ref(T... args)
get.hpp
Provides seqan3::views::get.
seqan3::field::alignment
@ alignment
The (pairwise) alignment stored in an seqan3::alignment object.
string