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;
}
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.";
}
};
Construction and assignment of alphabet symbols
using seqan3::operator""_dna4;
int main ()
{
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 it's 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>;
}
Reading records
After construction, you can now read the sequence records. Our file object behaves like a range so 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 & rec : 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
for (auto && records : fin | ranges::views::chunk(10))
{
}
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.
auto minimum_quality_filter = std::views::filter([] (auto const & rec)
{
auto qual = seqan3::get<seqan3::field::qual>(rec) |
std::views::transform([] (
auto q) {
return q.to_phred(); });
});
for (auto & rec : fin | minimum_quality_filter)
{
}
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.
for (auto && [rec1, rec2] : seqan3::views::zip(fin1, fin2))
{
if (seqan3::get<seqan3::field::id>(rec1) != seqan3::get<seqan3::field::id>(rec2))
}
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;
std::ranges::copy(fin, std::cpp20::back_inserter(records));
}
Note that you can move the record out of the file if you want to store it somewhere without copying.
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()
{
using seqan3::operator""_dna5;
for (int i = 0; i < 5; ++i)
{
seqan3::dna5_vector
seq{
"ACGT"_dna5};
fout.emplace_back(
seq,
id);
}
}
File conversion
Define a custom scoring scheme
using seqan3::operator""_dna4;
using seqan3::operator""_aa27;
auto sc_nc = nc_scheme.score('A'_dna4, 'C'_dna4);
auto sc_aa = aa_scheme.score('M'_aa27, 'K'_aa27);
- 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 less errors.
using seqan3::operator""_dna4;
int main()
{
"ACGAAGACCGAT"_dna4,
"ACGTGACTGACT"_dna4,
"AGGTACGAGCGACACT"_dna4};
auto filter_v = std::views::filter([](auto && res) { return res.score() >= -6;});
{
}
}
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.
using seqan3::operator""_dna4;
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();
}
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."};
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:
for (auto & [ id, alignment ] : fin)
{
}
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)
{
ref_seqs.push_back(
std::move(seqan3::get<seqan3::field::seq>(record)));
}
seqan3::field::ref_id,
#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 )
#else // ^^^ workaround / no workaround vvv
for (auto & [id, ref_id, mapq, alignment] : mapping_file | mapq_filter)
#endif // SEQAN3_WORKAROUND_GCC_93983
{
size_t sum_ref{};
for (auto const & char_ref : get<0>(alignment))
++sum_ref;
size_t sum_read{};
for (auto const & char_read : get<1>(alignment))
++sum_read;
<< sum_read << " gaps in the read sequence and "
<< sum_ref << " gaps in the reference sequence.\n";
}
}
Map reads ans 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);
}
seqan3::field::ref_id,
seqan3::field::ref_offset,
seqan3::field::qual,
for (auto && record : query_file_in)
{
auto & query = seqan3::get<seqan3::field::seq>(record);
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};
{
auto aligned_seq = alignment.alignment();
size_t ref_offset = alignment.sequence1_begin_position() + 2 + start;
size_t map_qual = 60u + alignment.score();
sam_out.emplace_back(query,
seqan3::get<seqan3::field::id>(record),
storage.ids[result.reference_id()],
ref_offset,
aligned_seq,
seqan3::get<seqan3::field::qual>(record),
map_qual);
}
}
}
}
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;
}
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)
{
bool flag{false};
top_level_parser.add_flag(flag, 'f', "flag", "some flag");
try
{
top_level_parser.parse();
}
{
return -1;
}
run_git_pull(sub_parser);
run_git_push(sub_parser);
else
}
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;
}
A custom dna4 alphabet that converts all unknown characters to <tt>A</tt>
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 fromchar` 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[
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;
}()
};
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);
}
{
'T'_my_dna4,
'G'_my_dna4,
'C'_my_dna4,
'A'_my_dna4
};
int main()
{
my_dna4 my_letter{'C'_my_dna4};
my_letter.assign_char('S');
}
If you are interested in custom alphabets, also take a look at our tutorial How to write your own alphabet.
![$\in[0,1]$](form_64.png)
1.8.20