SeqAn3 3.4.0-rc.1
The Modern C++ library for sequence analysis.
Loading...
Searching...
No Matches
Sequence File Input and Output

Learning Objective:
You will get an overview of how Input/Output files are handled in SeqAn and learn how to read and write sequence files. This tutorial is a walk-through with links into the API documentation and also meant as a source for copy-and-paste code.

DifficultyEasy
Duration90 min
Prerequisite tutorialsSetup, Ranges, Alphabets
Recommended reading

File I/O in SeqAn

Most file formats in bioinformatics are structured as lists of records. In SeqAn we model our files as a range over records. This interface allows us to easily stream over a file, apply filters and convert formats, sometimes only in a single line of code. The file format is automatically detected by the file name extension and compressed files can be handled without any effort. We can even stream over files in a python-like way with range-based for loops:

for (auto & record : file)
{
// do something with my record
}

We will explain the details about reading and writing files in the Sequence File section below. Currently, SeqAn supports the following file formats:

Warning
Access to compressed files relies on external libraries. For instance, you need to have zlib installed for reading .gz files and libbz2 for reading .bz2 files. You can check whether you have installed these libraries by running cmake . in your build directory. If -- Optional dependency: ZLIB-x.x.x found. is displayed on the command line then you can read/write compressed files in your programs.

Basic layout of SeqAn file objects

Before we dive into the details, we will outline the general design of our file objects, hoping that it will make the following tutorial easier to understand.

As mentioned above, our file object is a range over records. More specifically, a range over objects of type seqan3::sequence_record.

Output files can handle various types that fulfill the requirements of the format (e.g. a sequence has to be a range over an alphabet). In contrast to this, input files have certain default types for record fields that can be modified via a traits type.

For example, on construction you can specify seqan3::sequence_file_input_default_traits_dna or seqan3::sequence_file_input_default_traits_aa for reading dna (Nucleotide) and protein (Aminoacid) sequences respectively (section File traits will covers this in more detail).

Opening and closing files is also handled automatically. If a file cannot be opened for reading or writing, a seqan3::file_open_error is thrown.

Sequence file formats

Sequence files are the most generic and common biological files. Well-known formats include FASTA and FASTQ.

FASTA format

A FASTA record contains the sequence id and the sequence characters. Here is an example of a FASTA file:

>seq1
CCCCCCCCCCCCCCC
>seq2
CGATCGATC

In SeqAn we provide the seqan3::format_fasta to read sequence files in FASTA format.

FASTQ format

A FASTQ record contains an additional quality value for each sequence character. Here is an example of a FASTQ file:

@seq1
CCCCCCCCCCCCCCC
+
IIIIIHIIIIIIIII
@seq2
CGATCGATC
+
IIIIIIIII

In SeqAn we provide the seqan3::format_fastq to read sequence files in FASTQ format.

EMBL format

An EMBL record stores sequence and its annotation together. We are only interested in the id (ID) and sequence (SQ) information for a sequence file. Qualities are not stored in this format. Here is an example of an EMBL file:

ID X56734; SV 1; linear; mRNA; STD; PLN; 1859 BP.
XX
AC X56734; S46826;
XX
SQ Sequence 1859 BP; 609 A; 314 C; 355 G; 581 T; 0 other;
aaacaaacca aatatggatt ttattgtagc catatttgct ctgtttgtta ttagctcatt
cacaattact tccacaaatg cagttgaagc ttctactctt cttgacatag gtaacctgag

In SeqAn we provide the seqan3::format_embl to read sequence files in EMBL format.

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 SeqAn format class File Extensions
FASTA seqan3::format_fasta .fa, .fasta, .fna, .ffn, .ffa, .frn
FASTQ seqan3::format_fastq .fq, .fastq
EMBL seqan3::format_embl .embl

You can access the valid file extensions via the file_extensions member variable in a format and you can also customise this list if you want to allow different or additional file extensions:

int main()
{
seqan3::debug_stream << seqan3::format_fastq::file_extensions << '\n'; // prints [fastq,fq,qq]
seqan3::sequence_file_input fin{std::filesystem::current_path() / "my.qq"}; // detects FASTQ format
}
static std::vector< std::string > file_extensions
The valid file extensions for this format; note that you can modify this value.
Definition format_fastq.hpp:92
A class for reading sequence files, e.g. FASTA, FASTQ ...
Definition sequence_file/input.hpp:207
T current_path(T... args)
Provides seqan3::debug_stream and related types.
debug_stream_type debug_stream
A global instance of seqan3::debug_stream_type.
Definition debug_stream.hpp:37
Meta-header for the IO / Sequence File submodule .
T push_back(T... args)

Fields

The Sequence file abstraction supports reading three different fields:

  1. seqan3::sequence_record::sequence
  2. seqan3::sequence_record::id
  3. seqan3::sequence_record::base_qualities

The three fields are retrieved by default.

Reading a sequence file

You can include the SeqAn sequence file functionality with:

Construction

At first, you need to construct a seqan3::sequence_file_input object that handles file access. In most cases, you construct from a file name, but you can also construct a sequence file object directly from a stream (e.g. std::cin or std::stringstream), but then you need to know your format beforehand:

#include <filesystem>
int main()
{
auto filename = std::filesystem::current_path() / "my.fasta";
seqan3::sequence_file_input fin_from_filename{filename};
return 0;
}
The FASTA format.
Definition format_fasta.hpp:77

All template parameters of the seqan3::sequence_file_input are automatically deduced, even the format! We detect the format by the file name extension. The file extension in the example above is .fasta so we choose the seqan3::format_fasta.

File traits

The seqan3::sequence_file_input needs to know the types of the data you are reading in (e.g. that your sequence is a dna sequence) at compile-time. These necessary types are defined in the traits type argument (seqan3::sequence_file_input::traits_type) which by default is set to seqan3::sequence_file_input_default_traits_dna.

We thereby assume that

In case you want to read a protein sequence instead, we also provide the seqan3::sequence_file_input_default_traits_aa traits type which sets the seqan3::sequence_record::sequence field to a std::vector over the seqan3::aa27 alphabet.

You can specify the traits object as the first template argument for the sequence file:

int main()
{
}
A traits type that specifies input as amino acids.
Definition sequence_file/input.hpp:167

You can also customise the types by inheriting from one of the default traits and changing the type manually. See the detailed information on seqan3::sequence_file_input_default_traits_dna for an example.

Reading records

After construction, you can now read the seqan3::sequence_record's. As described in the basic file layout, our file objects behave like ranges so you can use a range based for loop to conveniently iterate over 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>
auto input = R"(>TEST1
ACGT
>Test2
AGGCTGA
>Test3
GGAGTATAATATATATATATATAT)";
int main()
{
for (auto & record : fin)
{
seqan3::debug_stream << "ID: " << record.id() << '\n';
seqan3::debug_stream << "SEQ: " << record.sequence() << '\n';
// a quality field also exists, but is not printed, because we know it's empty for FASTA files.
}
}
Provides seqan3::sequence_file_input and corresponding traits classes.
Attention
An input file is a single input range, which means you can only iterate over it once!

In the above example, record has the type seqan3::sequence_file_input::record_type which is seqan3::sequence_record.

Note
It is important to write auto & and not just auto, otherwise you will copy the record on every iteration.

You can read up more on the different ways to stream over the file object in the detailed documentation of seqan3::sequence_file_input.

Assignment 1: Reading a FASTQ file

Copy and paste the following FASTQ file to some location, e.g. the tmp directory:

@seq1
AGCTAGCAGCGATCG
+
IIIIIHIIIIIIIII
@seq2
CGATCGATC
+
IIIIIIIII
@seq3
AGCGATCGAGGAATATAT
+
IIIIHHGIIIIHHGIIIH

and then create a simple program that reads in all the records from that file and prints the id, sequence and quality information to the command line using the seqan3::debug_stream. Access the record via the member accessor of seqan3::sequence_record!

Note: you include the seqan3::debug_stream with the following header

Solution

#include <filesystem>
int main()
{
seqan3::sequence_file_input fin{current_path / "my.fastq"};
for (auto & rec : fin)
{
seqan3::debug_stream << "ID: " << rec.id() << '\n';
seqan3::debug_stream << "SEQ: " << rec.sequence() << '\n';
seqan3::debug_stream << "QUAL: " << rec.base_qualities() << '\n';
}
}

The code will print the following:

ID: seq1
SEQ: AGCTAGCAGCGATCG
QUAL: IIIIIHIIIIIIIII
ID: seq2
SEQ: CGATCGATC
QUAL: IIIIIIIII
ID: seq3
SEQ: AGCGATCGAGGAATATAT
QUAL: IIIIHHGIIIIHHGIIIH

The record type

In the examples above, we always use auto to deduce the record type automatically. In case you need the type explicitly, e.g. if you want to store the records in a variable, you can access the type member seqan3::sequence_file_input::record_type.

int main()
{
using record_type = typename decltype(fin)::record_type;
// Because `fin` is a range, we can access the first element by dereferencing fin.begin()
record_type rec = *fin.begin();
}
iterator begin()
Returns an iterator to current position in the file.
Definition sequence_file/input.hpp:411

You can move the record out of the file if you want to store it somewhere without copying.

int main()
{
using record_type = typename decltype(fin)::record_type;
record_type rec = std::move(*fin.begin()); // avoid copying
}

Assignment 2: Storing records in a std::vector

Create a small program that reads in a FASTA file and stores all the records in a std::vector.

After reading, print the vector (this works natively with the seqan3::debug_stream).

Test your program with the following file:

>seq1
AGCT
>seq2
CGATCGA

It should print the following:

[(AGCT,seq1,),(CGATCGA,seq2,)]

Note that the quality (third element in the triple) is empty because we are reading a FASTA file.

Solution

#include <filesystem>
#include <ranges> // std::ranges::copy
int main()
{
seqan3::sequence_file_input fin{current_path / "my.fasta"};
using record_type = decltype(fin)::record_type;
// You can use a for loop:
for (auto & record : fin)
{
records.push_back(std::move(record));
}
// But you can also do this:
seqan3::debug_stream << records << '\n';
}
T back_inserter(T... args)
T copy(T... args)

Sequence files as views

Since SeqAn files are ranges, you can also create views over files. This enables us to create solutions for lots of use cases using only a few lines of code.

Reading a file in chunks

A common use case is to read chunks from a file instead of the whole file at once or line by line.

You can do so easily on a file range by using the seqan3::views::chunk.

int main()
{
// `&&` is important because seqan3::views::chunk returns temporaries!
for (auto && records : fin | seqan3::views::chunk(10))
{
// `records` contains 10 elements (or less at the end)
seqan3::debug_stream << "Taking the next 10 sequences:\n";
seqan3::debug_stream << "ID: " << (*records.begin()).id() << '\n'; // prints first ID in batch
}
Provides seqan3::views::chunk.
The main SeqAn3 namespace.
Definition aligned_sequence_concept.hpp:26

The example above will iterate over the file by reading 10 records at a time. If no 10 records are available any more, it will just print the remaining records.

Applying a filter to a file

On some occasions, you are only interested in sequence records that fulfill a certain criteria, e.g. having a minimum sequence length or a minimum average quality. Just like in the example with seqan3::views::chunk you can use std::ranges::filter for this purpose:

#include <numeric> // std::accumulate
#include <ranges>
int main()
{
// std::views::filter takes a function object (a lambda in this case) as input that returns a boolean
auto minimum_quality_filter = std::views::filter(
[](auto const & rec)
{
auto qualities = rec.base_qualities()
| std::views::transform(
[](auto quality)
{
return seqan3::to_phred(quality);
});
auto sum = std::accumulate(qualities.begin(), qualities.end(), 0);
return sum / std::ranges::size(qualities) >= 40; // minimum average quality >= 40
});
for (auto & rec : fin | minimum_quality_filter)
{
seqan3::debug_stream << "ID: " << rec.id() << '\n';
}
}
T accumulate(T... args)
constexpr auto to_phred
The public getter function for the Phred representation of a quality score.
Definition alphabet/quality/concept.hpp:97

To remind you what you have learned in the Ranges tutorial before, a view is not applied immediately but lazy evaluated. That means your file is still parsed record by record and not at once.

Reading paired-end reads

In modern Next Generation Sequencing experiments, you often encounter paired-end read data which is split into two files. The read pairs are identified by their identical id and position in the files.

If you want to handle one pair of reads at a time, you can do so easily with a seqan3::views::zip.

int main()
{
// for simplicity we take the same file
for (auto && [rec1, rec2] : seqan3::views::zip(fin1, fin2)) // && is important!
{ // because seqan3::views::zip returns temporaries
if (rec1.id() != rec2.id())
throw std::runtime_error("Your pairs don't match.");
}
}
seqan::stl::views::zip zip
A view adaptor that takes several views and returns tuple-like values from every i-th element of each...
Definition zip.hpp:24
Provides seqan3::views::zip.

Assignment 3: Fun with file ranges

Implement a small program that reads in a FASTQ file and prints the first 2 sequences that have a length of at least 5.

Hints:

  • You can use std::ranges::size to retrieve the size of a range.
  • You need the following includes for std::views::filter and std::views::take
    #include <ranges>

Test your program on the following FASTQ file:

@seq1
CGATCGATC
+
IIIIIIIII
@seq2
AGCG
+
IIII
@seq3
AGCTAGCAGCGATCG
+
IIIIIHIIJJIIIII
@seq4
AGC
+
III
@seq5
AGCTAGCAGCGATCG
+
IIIIIHIIJJIIIII

It should output seq1 and seq3.

Solution

#include <filesystem>
#include <ranges>
int main()
{
seqan3::sequence_file_input fin{current_path / "my.fastq"};
auto length_filter = std::views::filter(
[](auto const & rec)
{
return std::ranges::size(rec.sequence()) >= 5;
});
// Store all IDs into a vector:
for (auto & record : fin | length_filter | std::views::take(2))
{
ids.push_back(std::move(record.id()));
}
seqan3::debug_stream << ids << '\n';
}
typename decltype(detail::split_after< i >(list_t{}))::first_type take
Return a seqan3::type_list of the first n types in the input type list.
Definition type_list/traits.hpp:374
SeqAn specific customisations in the standard namespace.

Writing a sequence file

Construction

You construct the seqan3::sequence_file_output just like the seqan3::sequence_file_input by giving it a file name or a stream.

#include <filesystem>
int main()
{
auto fasta_file = std::filesystem::current_path() / "my.fasta";
// FASTA format detected, std::ofstream opened for file
}
A class for writing sequence files, e.g. FASTA, FASTQ ...
Definition io/sequence_file/output.hpp:66
Provides seqan3::sequence_file_output and corresponding traits classes.

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
int main()
{
// ^ no need to specify the template arguments
}

Writing records

The easiest way to write to a sequence file is to use the seqan3::sequence_file_output::push_back() member function. It works similarly to how it works on a std::vector.

// 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 <string>
int main()
{
using namespace seqan3::literals;
using sequence_record_type = seqan3::sequence_record<types, fields>;
for (int i = 0; i < 5; ++i) // ...
{
std::string id{"test_id"};
seqan3::dna5_vector sequence{"ACGT"_dna5};
sequence_record_type record{std::move(sequence), std::move(id)};
fout.push_back(record);
}
}
The record type of seqan3::sequence_file_input.
Definition sequence_file/record.hpp:26
Provides seqan3::dna5, container aliases and string literals.
The generic concept for a (biological) sequence.
The SeqAn namespace for literals.
Provides seqan3::sequence_record.
A class template that holds a choice of seqan3::field.
Definition record.hpp:125
The class template that file records are based on; behaves like a std::tuple.
Definition record.hpp:190
Type that contains multiple types.
Definition type_list.hpp:26

Assignment 4: Writing a FASTQ file

Use your code (or the solution) from the previous assignment. Iterate over the records with a for loop and instead of just printing the ids, write out all the records that satisfy the filter to a new file called output.fastq.

Test your code on the same FASTQ file. The file output.fastq should contain the following records:

@seq1
CGATCGATC
+
IIIIIIIII
@seq2
AGCG
+
IIII
@seq3
AGCTAGCAGCGATCG
+
IIIIIHIIJJIIIII
@seq4
AGC
+
III
@seq5
AGCTAGCAGCGATCG
+
IIIIIHIIJJIIIII

Solution

#include <filesystem>
#include <ranges>
int main()
{
seqan3::sequence_file_input fin{current_path / "my.fastq"};
seqan3::sequence_file_output fout{current_path / "output.fastq"};
auto length_filter = std::views::filter(
[](auto const & rec)
{
return std::ranges::size(rec.sequence()) >= 5;
});
for (auto & record : fin | length_filter)
{
fout.push_back(record);
}
}

Files as views

Again we want to point out the convenient advantage of modelling files as ranges. The seqan3::sequence_file_input models std::ranges::input_range and in the "reading a file" section you already saw a few examples of how to pipe a view onto such a range. The output files model std::ranges::output_range. This allows us to use ranges algorithms like std::ranges::move to "move" records from input to output files.

Todo:
There is no standard way for piping output ranges (e.g. fin | fout;), yet. Reevaluate this section later.

In the same way you can pipe the output file:

#include <algorithm>
#include <ranges>
int main()
{
auto current_path = std::filesystem::current_path();
seqan3::sequence_file_input fin{current_path / "my.fastq"};
seqan3::sequence_file_output fout{current_path / "output.fastq"};
// the following are equivalent:
// 1. copy records of input file into output file
std::ranges::move(fin, fout.begin());
// 2. assign all records of input file to output file
fout = fin;
// 3. same as 2. but as one liner
}
T move(T... args)

Assignment 5: Fun with file ranges 2

Working on your solution from the previous assignment, try to remove the for loop in favour of a std::ranges algorithm.

The result should be the same.

Solution

#include <algorithm>
#include <filesystem>
#include <ranges>
int main()
{
seqan3::sequence_file_input fin{current_path / "my.fastq"};
seqan3::sequence_file_output fout{current_path / "output.fastq"};
auto length_filter = std::views::filter(
[](auto & rec)
{
return std::ranges::size(rec.sequence()) >= 5;
});
fout = fin | length_filter;
// This would also work: copy an input range into an output range
std::ranges::move(fin | length_filter, fout.begin());
}

File conversion

As mentioned before, the seqan3::sequence_record is shared between formats which allows for easy file conversion. For example you can read in a FASTQ file and output a FASTA file in one line:

int main()
{
auto current_path = std::filesystem::current_path();
seqan3::sequence_file_output{current_path / "output.fasta"} =
seqan3::sequence_file_input{current_path / "my.fastq"};
}

Yes, that's it! Of course this only works because all fields that are required in FASTA are provided in FASTQ. The other way around would not work as easily because we have no quality information (and would make less sense too).

Hide me