This document provides example recipes on how to carry out particular tasks using the SeqAn functionalities in C++. Please note that these recipes are not ordered. You can use the links in the table of contents or the search function of your browser to navigate them.
It will take some time, but we hope to expand this document into containing numerous great examples. If you have suggestions for how to improve the Cookbook and/or examples you would like included, please feel free to contact us.
Read sequence files
int main ()
{
for (auto & [seq, id, qual] : file_in)
{
}
return 0;
}
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
@ seq
The "sequence", usually a range of nucleotides or amino acids.
@ qual
The qualities, usually in Phred score notation.
T temp_directory_path(T... args)
Write a custom validator
This recipe implements a validator that checks if a numeric argument is an integral square (i.e. 0, 1, 4, 9...). Invalid values throw a seqan3::validation_error.
struct custom_validator
{
using option_value_type = double;
void operator() (option_value_type const & val) const
{
{
}
}
{
return "Value must be the square of an integral number.";
}
};
Argument parser exception thrown when an argument could not be casted to the according type.
Definition: exceptions.hpp:124
Construction and assignment of alphabet symbols
int main ()
{
Meta-header for the alphabet module.
constexpr derived_type & assign_char(char_type const chr) noexcept
Assign from a character, implicitly converts invalid characters.
Definition: alphabet_base.hpp:165
constexpr derived_type & assign_rank(rank_type const c) noexcept
Assign from a numeric value.
Definition: alphabet_base.hpp:191
The four letter DNA alphabet of A,C,G,T..
Definition: dna4.hpp:53
The SeqAn namespace for literals.
constexpr rank_type to_rank() const noexcept
Return the letter's numeric value (rank in the alphabet).
Definition: alphabet_base.hpp:139
decltype(seqan3::to_rank(std::declval< semi_alphabet_type >())) alphabet_rank_t
The rank_type of the semi-alphabet; defined as the return type of seqan3::to_rank....
Definition: concept.hpp:167
Reverse complement and the six-frame translation of a string using views
This recipe creates a small program that
- reads a string from the command line (first argument to the program)
- "converts" the string to a range of seqan3::dna5 (Bonus: throws an exception if loss of information occurs)
- prints the string and its reverse complement
- prints the six-frame translation of the string
int main(int argc, char** argv)
{
myparser.add_positional_option(s, "Please specify the DNA string.");
try
{
myparser.parse();
}
{
return 0;
}
auto s_as_dna = s | seqan3::views::char_to<seqan3::dna5>;
}
Meta-header for the Alphabet / Views submodule .
Meta-header for the Argument Parser module .
Argument parser exception that is thrown whenever there is an error while parsing the command line ar...
Definition: exceptions.hpp:40
The SeqAn command line parser.
Definition: argument_parser.hpp:154
constexpr auto translate
A view that translates nucleotide into aminoacid alphabet with 1, 2, 3 or 6 frames.
Definition: translate.hpp:813
auto const complement
A view that converts a range of nucleotides to their complement.
Definition: complement.hpp:67
The <ranges> header from C++20's standard library.
Reading records
After construction, you can now read the sequence records. Our file object behaves like a range, you can use a range-based for loop to conveniently iterate over the file:
auto input = R"(> TEST1
ACGT
> Test2
AGGCTGA
> Test3
GGAGTATAATATATATATATATAT)";
int main()
{
for (auto & record : fin)
{
}
}
- Attention
- An input file is a single input range, which means you can only iterate over it once!
- Note
- It is important to write
auto &
and not just auto
, otherwise you will copy the record on every iteration.
You can also use structured binding, i.e. for (auto & [seq, id, qual] : fin)
But beware: with structured bindings you do need to get the order of elements correct!
You can also read a file in chunks:
Reading records in chunks
int main()
{
{
}
Provides seqan3::views::chunk.
T current_path(T... args)
constexpr auto chunk
A chunk view.
Definition: chunk.hpp:29
Meta-header for the IO / Sequence File submodule .
The example above will iterate over the file by reading 10 records at a time. If no 10 records are available anymore, 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 criterion, e.g. having a minimum sequence length or a minimum average quality.
This recipe can be used to filter the sequences in your file by a custom criterion.
int main()
{
auto minimum_quality_filter = std::views::filter([] (auto const & rec)
{
auto qualities = rec.base_qualities()
});
for (auto & rec : fin | minimum_quality_filter)
{
}
}
constexpr auto to_phred
The public getter function for the Phred representation of a quality score.
Definition: concept.hpp:100
decltype(detail::transform< trait_t >(list_t{})) transform
Apply a transformation trait to every type in the list and return a seqan3::type_list of the results.
Definition: traits.hpp:471
constexpr size_t size
The size of a type pack.
Definition: traits.hpp:151
Reading paired-end reads
In modern Next Generation Sequencing experiments you often have paired-end read data which is split into two files. The read pairs are identified by their identical name/id and position in the two files.
This recipe can be used to handle one pair of reads at a time.
int main()
{
{
if (rec1.id() != rec2.id())
}
}
constexpr auto zip
A zip view.
Definition: zip.hpp:29
Provides seqan3::views::zip.
Storing records in a std::vector
This recipe creates a small program that reads in a FASTA file and stores all the records in a std::vector.
int main()
{
using record_type = decltype(fin)::record_type;
for (auto & record : fin)
{
}
std::ranges::copy(fin, std::cpp20::back_inserter(records));
}
The <filesystem> header from C++17's standard library.
Note that 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());
}
Writing records
The easiest way to write to a sequence file is to use the seqan3::sequence_file_output::push_back() or seqan3::sequence_file_output::emplace_back() member functions. These work similarly to how they work on a std::vector.
int main()
{
for (int i = 0; i < 5; ++i)
{
seqan3::dna5_vector
sequence{
"ACGT"_dna5};
sequence_record_type record{std::move(
sequence), std::move(
id)};
fout.push_back(record);
}
}
A class for writing sequence files, e.g. FASTA, FASTQ ...
Definition: output.hpp:69
The record type of seqan3::sequence_file_input.
Definition: record.hpp:29
Provides seqan3::dna5, container aliases and string literals.
The generic concept for a (biological) sequence.
Provides seqan3::sequence_file_output and corresponding traits classes.
Provides seqan3::sequence_record.
A class template that holds a choice of seqan3::field.
Definition: record.hpp:128
Type that contains multiple types.
Definition: type_list.hpp:29
File conversion
Define a custom scoring scheme
Provides seqan3::aminoacid_scoring_scheme.
Provides seqan3::nucleotide_scoring_scheme.
auto sc_nc = nc_scheme.score('A'_dna4, 'C'_dna4);
auto sc_aa = aa_scheme.score('M'_aa27, 'K'_aa27);
A data structure for managing and computing the score of two amino acids.
Definition: aminoacid_scoring_scheme.hpp:75
constexpr void set_similarity_matrix(aminoacid_similarity_matrix const matrix_id)
Set the similarity matrix scheme (e.g. blosum62).
Definition: aminoacid_scoring_scheme.hpp:122
A data structure for managing and computing the score of two nucleotides.
Definition: nucleotide_scoring_scheme.hpp:38
@ blosum30
The blosum30 matrix for very distantly related proteins.
A strong type of underlying type score_type that represents the score of two matching characters.
Definition: scoring_scheme_base.hpp:41
A strong type of underlying type score_type that represents the score two different characters.
Definition: scoring_scheme_base.hpp:66
- Attention
- SeqAn's alignment algorithm computes the maximal similarity score, thus the match score must be set to a positive value and the scores for mismatch and gap must be negative in order to maximize over the matching letters.
Calculate edit distance for a set of sequences
This recipe can be used to calculate the edit distance for all six pairwise combinations. Here we only allow at most 7 errors and filter all alignments with 6 or fewer errors.
int main()
{
"ACGAAGACCGAT"_dna4,
"ACGTGACTGACT"_dna4,
"AGGTACGAGCGACACT"_dna4};
auto filter_v = std::views::filter([](auto && res) { return res.score() >= -6;});
for (auto const & result : alignment_results | filter_v)
{
}
}
Provides pairwise alignment function.
Sets the global alignment method.
Definition: align_config_method.hpp:122
Sets the minimal score (maximal errors) allowed during an distance computation e.g....
Definition: align_config_min_score.hpp:39
Configures the alignment result to output the score.
Definition: align_config_output.hpp:43
Provides seqan3::dna4, container aliases and string literals.
constexpr configuration edit_scheme
Shortcut for edit distance configuration.
Definition: align_config_edit.hpp:51
constexpr auto align_pairwise(sequence_t &&seq, alignment_config_t const &config)
Computes the pairwise alignment for a pair of sequences or a range over sequence pairs.
Definition: align_pairwise.hpp:135
constexpr auto pairwise_combine
A view adaptor that generates all pairwise combinations of the elements of the underlying range.
Definition: pairwise_combine.hpp:710
Provides seqan3::views::pairwise_combine.
Searching for matches
This recipe can be used to search for all occurrences of a substring and print the number of hits and the positions in an ascending ordering.
void run_text_single()
{
seqan3::dna4_vector
text{"CGCTGTCTGAAGGATGAGTGTCAGCCAGTGTAACCCGATGAGCTACCCAGTAGTCGAACTGGGCCAGACAACCCGGCGCTAATGCACTCA"_dna4};
<< "The following hits were found:\n";
for (
auto && result :
search(
"GCT"_dna4, index))
}
void run_text_collection()
{
"ACCCGATGAGCTACCCAGTAGTCGAACTG"_dna4,
"GGCCAGACAACCCGGCGCTAATGCACTCA"_dna4};
<< "The following hits were found:\n";
for (
auto && result :
search(
"GCT"_dna4, index))
}
int main()
{
run_text_single();
run_text_collection();
}
The SeqAn FM Index.
Definition: fm_index.hpp:192
Provides the unidirectional seqan3::fm_index.
auto search(queries_t &&queries, index_t const &index, configuration_t const &cfg=search_cfg::default_configuration)
Search a query or a range of queries in an index.
Definition: search.hpp:104
Provides the public interface for search algorithms.
If you want to allow errors in your query, you need to configure the approximate search with the following search configuration objects:
To search for either 1 insertion or 1 deletion you can use the seqan3::search_cfg::error_count:
std::string text{
"Garfield the fat cat without a hat."};
Collection of elements to configure an algorithm.
Definition: configuration.hpp:45
Configuration element that represents the number or rate of deletion errors.
Definition: max_error.hpp:173
Configuration element that represents the number or rate of insertion errors.
Definition: max_error.hpp:127
Configuration element that represents the number or rate of substitution errors.
Definition: max_error.hpp:82
Configuration element that represents the number or rate of total errors.
Definition: max_error.hpp:37
A strong type of underlying type uint8_t that represents the number of errors.
Definition: max_error_common.hpp:32
Reading the CIGAR information into an actual alignment
In SeqAn, the conversion from a CIGAR string to an alignment (two seqan3::aligned_sequence
s) is done automatically for you. You can access the resulting alignment via seqan3::sam_record::alignment:
int main()
{
for (auto & record : fin)
}
Meta-header for the IO / SAM File submodule .
Combining sequence and alignment files
This recipe can be used to:
- Read in a FASTA file with the reference and a SAM file with the alignment
- Filter the alignment records and only take those with a mapping quality >= 30.
- For the resulting alignments, print which read was mapped against with reference id and the number of seqan3::gap's involved in the alignment (either in aligned reference or in read sequence).
int main()
{
for (auto && record : reference_file)
{
reference_ids.
push_back(std::move(record.id()));
reference_sequences.push_back(std::move(record.sequence()));
}
auto mapq_filter = std::views::filter([] (auto & record) { return record.mapping_quality() >= 30; });
for (auto & record : mapping_file | mapq_filter)
{
size_t sum_reference{};
for (auto const & char_reference : std::get<0>(record.alignment()))
++sum_reference;
<< " with "
<< sum_read << " gaps in the read sequence and "
<< sum_reference << " gaps in the reference sequence.\n";
}
}
The <algorithm> header from C++20's standard library.
The alphabet of a gap character '-'.
Definition: gap.hpp:39
constexpr ptrdiff_t count
Count the occurrences of a type in a pack.
Definition: traits.hpp:169
Provides the seqan3::record template and the seqan3::field enum.
Map reads and write output to SAM file
For a full recipe on creating your own readmapper, see the very end of the tutorial Implementing your own read mapper with SeqAn.
reference_storage_t & storage,
uint8_t const errors)
{
{
cereal::BinaryInputArchive iarchive{is};
iarchive(index);
}
for (auto && record : query_file_in)
{
auto & query = record.sequence();
for (
auto && result :
search(query, index, search_config))
{
size_t start = result.reference_begin_position() ? result.reference_begin_position() - 1 : 0;
std::span text_view{
std::data(storage.seqs[result.reference_id()]) + start, query.size() + 1};
{
sam_out.emplace_back(query,
record.id(),
storage.ids[result.reference_id()],
ref_offset,
aligned_seq,
record.base_qualities(),
map_qual);
}
}
}
}
Configures the alignment result to output the alignment.
Definition: align_config_output.hpp:171
Configures the alignment result to output the begin positions.
Definition: align_config_output.hpp:131
The SeqAn Bidirectional FM Index.
Definition: bi_fm_index.hpp:61
A class for writing alignment files, e.g. SAM, BAL, BLAST, ...
Definition: output.hpp:77
Configuration element to receive all hits with the lowest number of errors within the error bounds.
Definition: hit.hpp:59
@ ref_offset
Sequence (seqan3::field::ref_seq) relative start position (0-based), unsigned value.
@ alignment
The (pairwise) alignment stored in an object that models seqan3::detail::pairwise_alignment.
@ mapq
The mapping quality of the seqan3::field::seq alignment, usually a Phred-scaled score.
@ ref_id
The identifier of the (reference) sequence that seqan3::field::seq was aligned to.
@ id
The identifier, usually a string.
A strong type representing free_end_gaps_sequence1_leading of the seqan3::align_cfg::method_global.
Definition: align_config_method.hpp:68
A strong type representing free_end_gaps_sequence1_trailing of the seqan3::align_cfg::method_global.
Definition: align_config_method.hpp:88
A strong type representing free_end_gaps_sequence2_leading of the seqan3::align_cfg::method_global.
Definition: align_config_method.hpp:78
A strong type representing free_end_gaps_sequence2_trailing of the seqan3::align_cfg::method_global.
Definition: align_config_method.hpp:98
Constructing a basic argument parser
{
}
struct cmd_arguments
{
};
{
parser.
add_option(args.reference_path,
'r',
"reference",
"The path to the reference.",
parser.
add_option(args.index_path,
'o',
"output",
"The output index file path.",
}
int main(int argc, char const ** argv)
{
cmd_arguments args{};
initialise_argument_parser(parser, args);
try
{
parser.parse();
}
{
return -1;
}
run_program(args.reference_path, args.index_path);
return 0;
}
argument_parser_meta_data info
Aggregates all parser related meta data (see seqan3::argument_parser_meta_data struct).
Definition: argument_parser.hpp:604
void add_option(option_type &value, char const short_id, std::string const &long_id, std::string const &desc, option_spec const spec=option_spec::standard, validator_type option_validator=validator_type{})
Adds an option to the seqan3::argument_parser.
Definition: argument_parser.hpp:247
A validator that checks if a given path is a valid output file.
Definition: validators.hpp:646
@ standard
The default were no checking or special displaying is happening.
Definition: auxiliary.hpp:250
@ required
Definition: auxiliary.hpp:251
Constructing a subcommand argument parser
struct pull_arguments
{
bool progress{false};
};
{
pull_arguments args{};
try
{
}
{
return -1;
}
seqan3::debug_stream <<
"Git pull with repository " << args.repository <<
" and branch " << args.branch <<
'\n';
return 0;
}
struct push_arguments
{
bool push_all{false};
};
{
push_arguments args{};
parser.
add_positional_option(args.branches,
"The branch names to push (if none are given, push current).");
try
{
}
{
return -1;
}
seqan3::debug_stream <<
"Git push with repository " << args.repository <<
" and branches " << args.branches <<
'\n';
return 0;
}
int main(int argc, char const ** argv)
{
argc,
argv,
{"push", "pull"}};
top_level_parser.info.description.push_back("You can push or pull from a remote repository.");
top_level_parser.add_flag(flag, 'f', "flag", "some flag");
try
{
top_level_parser.parse();
}
{
return -1;
}
return run_git_pull(sub_parser);
return run_git_push(sub_parser);
else
return 0;
}
void add_positional_option(option_type &value, std::string const &desc, validator_type option_validator=validator_type{})
Adds a positional option to the seqan3::argument_parser.
Definition: argument_parser.hpp:310
void parse()
Initiates the actual command line parsing.
Definition: argument_parser.hpp:395
argument_parser & get_sub_parser()
Returns a reference to the sub-parser instance if subcommand parsing was enabled.
Definition: argument_parser.hpp:423
@ flag
The alignment flag (bit information), uint16_t value.
@ on
Automatic update notifications should be enabled.
Serialise data structures with cereal
#include <cereal/archives/binary.hpp>
#include <cereal/types/vector.hpp>
#include <seqan3/test/tmp_filename.hpp>
{
cereal::BinaryInputArchive archive(is);
archive(data);
}
{
cereal::BinaryOutputArchive archive(os);
archive(data);
}
int main()
{
seqan3::test::tmp_filename tmp_file{"data.out"};
store(vec, tmp_file);
load(vec2, tmp_file);
return 0;
}
Converting a range of an alphabet
using seqan3::operator""_dna4;
using seqan3::operator""_dna5;
using seqan3::operator""_phred42;
int main()
{
{'C'_dna4, 'A'_phred42},
{'G'_dna4, '6'_phred42},
{'T'_dna4, '&'_phred42}};
auto view2 = sequence2
{
return static_cast<seqan3::dna4>(std::forward<decltype(in)>(in));
});
return 0;
}
Provides aliases for qualified.
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: traits.hpp:368
Provides seqan3::phred42 quality scores.
Provides quality alphabet composites.
A custom dna4 alphabet that converts all unknown characters to A
When assigning from char
or converting from a larger nucleotide alphabet to a smaller one, loss of information can occur since obviously some bases are not available. When converting to seqan3::dna5 or seqan3::rna5, non-canonical bases (letters other than A, C, G, T, U) are converted to 'N'
to preserve ambiguity at that position. For seqan3::dna4 and seqan3::rna4 there is no letter 'N'
to represent ambiguity, so the conversion from char
for IUPAC characters tries to choose the best fitting alternative (see seqan3::dna4 for more details).
If you would like to always convert unknown characters to A
instead, you can create your own alphabet with a respective char conversion table very easily like this:
{
public:
using nucleotide_base<my_dna4, 4>::nucleotide_base;
private:
static constexpr char_type rank_to_char(rank_type const rank)
{
return rank_to_char_table[rank];
}
static constexpr rank_type char_to_rank(char_type const chr)
{
return char_to_rank_table[static_cast<index_t>(chr)];
}
static constexpr rank_type rank_complement(rank_type const rank)
{
return rank_complement_table[rank];
}
private:
static constexpr char_type rank_to_char_table[
alphabet_size] {
'A',
'C',
'G',
'T'};
{
[] () constexpr
{
conversion_table['C'] = conversion_table['c'] = 1;
conversion_table['G'] = conversion_table['g'] = 2;
conversion_table['T'] = conversion_table['t'] = 3;
conversion_table['U'] = conversion_table['T'];
conversion_table['u'] = conversion_table['t'];
return conversion_table;
}()
};
{
3,
2,
1,
0
};
friend nucleotide_base<my_dna4, 4>;
friend nucleotide_base<my_dna4, 4>::base_t;
};
constexpr my_dna4 operator""_my_dna4(char const c) noexcept
{
return my_dna4{}.assign_char(c);
}
int main()
{
my_dna4 my_letter{'C'_my_dna4};
my_letter.assign_char('S');
}
static constexpr detail::min_viable_uint_t< size > alphabet_size
The size of the alphabet, i.e. the number of different values it can take.
Definition: alphabet_base.hpp:203
A CRTP-base that refines seqan3::alphabet_base and is used by the nucleotides.
Definition: nucleotide_base.hpp:43
constexpr auto complement
Return the complement of a nucleotide object.
Definition: concept.hpp:104
Provides seqan3::nucleotide_base.
If you are interested in custom alphabets, also take a look at our tutorial How to write your own alphabet.
Controlling threads of (de-)compression streams
When reading or writing compressed files, parallelisation is automatically applied when using BGZF-compressed files, e.g., BAM files. This will use all available threads by default and can be adjusted by setting seqan3::contrib::bgzf_thread_count
to the desired value:
int main()
{
seqan3::contrib::bgzf_thread_count = 1u;
return 0;
}
Meta-header for the IO module .