Well, as SeqAn is written in C++ it might not come as a surprise, that you should bring some basic knowledge about the C++ programming language.
We made quite an effort to not depend on other libraries, than the on-board tools already bring.
This means, to learn SeqAn it suffices in the beginning to only know about C++ and the STL which is the standard template library defined by the ISO C++ committee.
The rest will be discussed in the subsequent tutorials step by step.

Before we start, here is a strong advice!
If you are diving into C++ for the first time, because you are new to programming or switched from another programming language, then we recommend, you first stroll through the C++ FAQs to acquaint yourself with C++.
There you can find many useful tips about C++ and get some further readings.
It also will introduce you to the paradigms that we used for designing this library.
In the Getting Started section we will introduce you to the design decisions of SeqAn and lay down some of the basic programming paradigms we follow to make this library so efficient.
If you never heard about these paradigms, dont‚Äôt worry.
We will give you code examples, which you can try out on your own.
Never forgot, there is no better way to learn a new language or language feature, than to actually program with it.
So keep your fingers attached to the keyboard and let‚Äôs start right away!

If you didn‚Äôt install SeqAn yet, please follow the User Guide instructions to install SeqAn first.
After that you should continue with the tutorials.

Hint

Please note, that although we try hard to provide a very comprehensive list of topics, it is not always possible to cover every angle of the library and its features.
The tutorials are thought as a first place to start.
If you are more experienced with SeqAn you can use the API documentation in addition to search for specific functions or classes.

The tutorial section is organized such that you can efficiently search for a specific topic you want to learn more about.
Each tutorial takes 30 to 60 minutes of your time for learning how to use SeqAn.
So buckle up and jump right into using SeqAn using our tutorials!

These articles are required for every one that is new to SeqAn.
Take your time and study these documents thoroughly, as they describe the fundamental concepts and design decisions of the library.
Everything else depends on these informations.

In this section we explain several different algorithms that are crucial for many bioinformatics applications.
This includes pattern matching, dynamic programming algorithms for sequence alignments, seed extension and many more.
Beginners that come from the tutorials about data structures should either continue with Online Pattern Matching or with the DP Alignment Algorithms.

On this page you will learn how to read/write and work with common bioinformatic file formats, such as FASTA, BAM, BED, VCF files, and more.
Beginners should start with the File I/O Overview.
This tutorial introduces you to the basic I/O concepts and data structures.

The how-to page is divided into Recipes and Use Cases.
The former section gives you some useful hints about miscellaneous topics.
The latter section describes how some use cases can be solved with SeqAn.
Things presented here are for experienced SeqAn users.
If you are a beginner, first have a look at the tutorials above.

You will learn about the design goals and fundamental ideas used in the SeqAn library.
Also, you will see how the SeqAn library can be generic while still retaining high performance.

Difficulty

Very basic

Duration

Take the time you need!

Prerequisites

Basic C or C++ knowledge

Hi, we are glad you made it here.
You being here, and reading these lines means you are eager to learn more about SeqAn and this is the right place to start.
In this tutorial, we will give you an overview about the design goals, design decisions of the SeqAn library, and explain the motivation for these decisions.
The next chapter First Steps will flesh out the most important points of this chapter with code examples of everyday SeqAn use.

The following lists some library design aims of the SeqAn library.
Note that they are contradicting.
The focus is on efficiency but small trade-offs are allowed to improve consistency and ease of use.

Efficiency.
The focus of SeqAn is to provide a library of efficient and reusable algorithmic components for biological sequence analysis.
Algorithms should have good practical implementations with low overhead, even at the cost of being harder to use.

Consistency.
Be consistent wherever possible, even at slight costs of efficiency.

Ease of use.
The library should be easily usable wherever possible, even at slight costs of efficiency.

Reusability and Generosity.
The algorithms in SeqAn should be reusable and generic, even at small costs of efficiency.

C++ is sometimes described as a language that most people know only 20% of but everyone knows a different 20%.
This section gives an overview over some C++ idioms we use.
This might be no news if you are a seasoned C++ programmer who is apt at using the STL and Boost libraries.
However, programmers coming from C and Java might find them interesting (We still encourage to read the C++ FAQ if you are new to C++).

C++ allows you to perform generic programming using templates.
While similar to generics in Java (C++ templates are more than a decade older), C++ templates are designed to write zero-overhead abstractions that can be written to be as efficient as hand-written code while retaining a high level of abstraction.
See cplusplus.com on C++ templates.
Note that there is no way to restrict the type that can be used in templates, there is no mechanism such as Java‚Äôs ?extendsT in C++.
Using an incompatible type leads to compiler errors because some operator or function could not be found.

Memory Management / No Pointers

Object oriented programming is another key programming paradigm made available with C++ (Compared to C).
This means, that instead of using raw pointers to allocated chunks of memory, memory management should be done using containers.
The STL provides containers such as std::vector and SeqAn offers String.

C++ allows to allocate complex objects on the stack (in contrast to Java where objects are always constructed on the heap).
The objects are constructed when the code execution enters the scope/block they are defined in and freed when the block is left.
Allocation of resources (e.g. memory) happens on construction and deallocation happens when the current block is left.
This is best explained in an example.

seqan::String<char> is a class (actually an instantiation of the class template String) that allows to store strings of char values, similar to std::vector<char> or std::string.

When the variable programName is allocated, the constructor of the String<char> class is called.
It allocates sufficient memory to store the value of argv[0] and then copies over the values from this string.
The variable exists until the current block is left.
Since it is defined in the main() function, this can only happen in the last line of main() at the return0.
When the variable goes out of scope, its value is deconstructed and all allocated memory is freed.

If an argument was given to the program, the block in the if clause is entered.
When this happens, the variable firstArg is constructed, memory is allocated and the value of argv[1] is copied into the buffer.
When the block is left, the variable is deconstructed and all memory is deallocated.

Note that all memory is released when the main() function is left, regardless whether it is left in the return0 or the return1.
Corresponding code in C would be (arguably) more messy, either requiring goto or multiple free() calls, one before either return.

In this section, we will give a short rationale why C++ with heavy use of template programming was used for SeqAn.

Any sequence analysis will have sequence data structures and algorithms on sequences at its heart.
Even when only considering DNA and amino acid alphabets, there are various variants for alphabets that one has to consider.
Otherwise, important applications in bioinformatics cannot be covered:

4-character DNA,

5-character DNA with N,

15-character IUPAC, and

27-character amino acids.

A simple implementation could simply store such strings as ASCII characters.
However, there are some implementation tricks that can lead to great reduction of memory usage (e.g. encoding eight 4-character DNA characters in one byte) or running time (fast lookup tables for characters or q-grams) for small alphabets.
Thus, simply using a std::string would come at high costs to efficiency.

Given that in the last 10-15 years, Java and C# have gained popularity, one could think about an object oriented solution: strings could simply be arrays of Character objects.
Using polymorphism (e.g. overwriting of functions in subclasses), one could then write generic and reusable algorithms.
For example, the Java 2 platform defines the sort function for all objects implementing a Comparable interface.
Note that such an implementation would have to rely on virtual functions of some sort.
However, as we will see in the section OOP vs. Generic Progamming, this comes at a high performance cost, being in conflict with efficiency.
For a sequence library, we could implement functions that map values from an alphabet to an ordinal value between 0 and S-1 where S is the number of elements in the alphabet.

Generic programming offers one way out: C++ templates allow to define template classes, e.g. the STL‚Äôs std::vector<T> or SeqAn‚Äôs String.
Here, instead of creating a string class around an array of char values (or objects), we can leave the type of the array‚Äôs elements open.
We can then introduce different types, e.g. Dna or Dna5 for 4- and 5-character DNA alphabets.

Algorithms can be implemented using template functions and the template types are fixed at compile time.
Thus, the compiler does not have to use virtual function tables and other ‚Äúcrutches‚ÄĚ, less indirection is involved, and more code can be inlined and aggressively optimized.
When written appropriately, such algorithms can also work on different string implementations! Also, when defining our own alphabet types, we can directly influence how their abstractions (and APIs) work.

Thus, C++ allows us to implement (1) a generic and reusable library with (2) high level abstractions (and thus ease of use) that still allows the compiler to employ aggressive optimization and thus achieves (3) efficiency.
With the words of the C++ inventor Bjarne Stroustrup:

A high level of abstraction is good, not just in C++, but in general.
We want to deal with problems at the level we are thinking about those problems.
When we do that, we have no gap between the way we understand problems and the way we implement their solutions.
We can understand the next guy‚Äôs code. We don‚Äôt have to be the compiler.

In SeqAn, we use a technique called template subclassing which is based on generic programming.
This technique provides polymorphism into C++ programs at compile time using templates.
Such static polymorphism is different from runtime polymorphism which is supported in C++ using subclassing and virtual functions.
It comes at the cost of some additional typing but has the advantage that the compiler can inline all function calls and thus achieve better performance.
An example will be given in the section ‚ÄúFrom OOP to SeqAn‚ÄĚ in the First Example Tutorial.

Todo

We need a little code example here.

The important point is that in contrast to runtime polymorphism such static polymorphism allows the compiler to inline functions, which has huge effect on the overall performance of the program.
Which as you recall correctly from above, is the main objective of the SeqAn library :)!

As we already stated, using template subclassing to achieve OOP like behavior in a more efficient way comes with a certain drawback.
Subclassed objects are seen by the compiler as singular instances of a specific type.
That means a subclassed object does not inherit the member or member functions of the alleged base class.
In order to reduce the overhead of reimplementing the same member functions for every subclassed object, we use global interface functions.

You might already have seen global function interfaces while working with the STL.
With the new C++11 standard the STL now provides some global interface functions, e.g., the begin or end interface.

The rationale behind this is the following observation.
Global interface functions allow us to implement a general functionality that is used for all subclassed objects of this template class (assuming the accessed member variables exists in all subclassed objects as in the base template class, otherwise the compiler will complain).
If the behavior for any subclassed object changes, the corresponding global function will be reimplemented for this special type covering the desired functionality.
Due to template deduction the compiler already chooses the correct function and inlines the kernel if possible, which very likely improves the performance of the program.
By this design, we can avoid code duplication, and by that increasing maintainability and reducing subtle errors due to less copy-and-paste code.

So, while most C++ developers, who are familiar with the STL and have a strong background in OO programming, are used to the typical dot notation, in SeqAn you have to get used to global function interfaces instead.
But, cheer up! You will adapt to this very quickly. Promised!

Generic algorithms usually have to know certain types that correspond to their arguments.
An algorithm on containers may need to know which type of values are stored in the string, or what kind of iterator we need to access it.
The usual way in the STL is to define the value type of a class like vector as a member typedef of this class, so it can be retrieved by vector::value_type.

Unfortunately member typedef declarations have the same disadvantages as any members: Since they are specified by the class definition, they cannot be changed or added to the class without changing the code of the class, and it is not possible in C++ to define members for built-in types.
What we need therefore is a mechanism that returns an output type (e.g. the value type) given an input type (e.g. the string) and doing so does not rely on members of the input type, but instead uses some kind of global interface.

Such tasks can be performed by metafunctions, also known as type traits.
A metafunction is a construct to map some types or constants to other entities like types, constants, functions, or objects at compile time.
The name metafunction comes from fact that they can be regarded as part of a meta-programming language that is evaluated during compilation.

In SeqAn we use class templates to implement metafunctions in C++.
Generic algorithms usually have to know certain types that correspond to their arguments: An algorithm on strings may need to know which type of characters are stored in the string, or what kind of iterator can be used to browse it.
SeqAn uses Metafunctions (also known as ‚Äútraits‚ÄĚ) for that purpose.

Now the function works for all sequence types T that store AminoAcid objects, but it will fail for other value types as soon as the variable temp cannot store str[0] anymore.
To overcome this problem, we must redefine temp in a way that it can store a value of the correct type.
The question is: ‚ÄúGiven a arbitrary type T, what is the value type of T?‚ÄĚ

The metafunction Value answers this question: ‚ÄúThe value type of T is given by Value<T>::Type.‚ÄĚ

Hence, the final version of our function exchangeFirstValues reads as follows:

We can view Value as a kind of ‚Äúfunction‚ÄĚ that takes T as an argument (in angle brackets) and returns the required value type of T.
In fact, Value is not implemented as a C++ function, but as a class template.
This class template is specialized for each sequence type T in a way that the typedefType provides the value type of T.
Unfortunately, the current C++ language standard does not allow to write simply ‚ÄúValue<T>temp;‚ÄĚ, so we must select the return value by appending ‚Äú::Type‚ÄĚ.
The leading ‚Äútypename‚ÄĚ becomes necessary since Value<T>::Type is a type that depends on a template parameter of the surrounding function template.

Wow, this was quite some information to digest, wasn‚Äôt it?
We suggest you take a break!
Get some fresh air!
Grab something to drink or to eat!
Let the information settle down.

Do you think you‚Äôve got everything?
Well, if not don‚Äôt worry!
Follow the First Steps tutorial which will cover the topics discussed above.
This gives you the chance to apply the recently discussed paradigms to an actual (uhm, simplistic) use case.
But it will help you to better understand the way data structures and algorithms are implemented in SeqAn.

We recommend you to also read the Argument Parser Tutorial.
This tutorial will teach you how to easily add command line arguments for your program and how to generate a help page for the options.
Or you go back to the main page and stroll through the other tutorials.
You are now ready to dive deeper into SeqAn.
Enjoy!

You will learn the most basic concepts of SeqAn.
After this tutorial you will be ready to deal with the more specific tutorials, e.g. Sequences.

Difficulty

Very basic

Duration

1.5h

Prerequisites

Basic C or C++ knowledge

Welcome to the SeqAn ‚ÄúHello World‚ÄĚ.
This is the first practical tutorial you should look at when starting to use our software library.

We assume that you have some programming experience (preferably in C++ or C) and concentrate on SeqAn specific aspects.
We will start out pretty slowly and hopefully the tutorial will make sense to you even if you are new to C++.
However, to really leverage the power of SeqAn you will have to learn C++.
There are many tutorials on C++, for example the tutorial at cplusplus.com.

This tutorial will walk you through a simple example program that highlights the things that are most prominently different from the libraries that many SeqAn newcomers are used to:

Let‚Äôs start with a simple example programm. The program will do a pattern search of a short query sequence (pattern) in a long subject sequence (text).
We define the score for each position of the database sequence as the sum of matching characters between the pattern and the text.

The first position has a score of 1, because the i in the pattern matches the i in is.
This is only a toy example for explanatory reasons and we ignore any more advanced implementations.

In SeqAn the program could look like this (we will explain every line of code shortly):

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>usingnamespaceseqan;intmain(){// InitializationString<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score;resize(score,length(text)-length(pattern)+1);// Computation of the similarities// Iteration over the text (outer loop)for(unsignedi=0;i<length(text)-length(pattern)+1;++i){intlocalScore=0;// Iteration over the pattern for character comparisonfor(unsignedj=0;j<length(pattern);++j){if(text[i+j]==pattern[j])++localScore;}score[i]=localScore;}// Printing the resultfor(unsignedi=0;i<length(score);++i)std::cout<<score[i]<<" ";std::cout<<std::endl;return0;}

Whenever we use SeqAn classes or functions we have to explicitly write the namespace qualifier seqan:: in front of the class name or function.
This can be circumvented if we include the line usingnamespaceseqan; at the top of the working example.
However, during this tutorial we will not do this, such that SeqAn classes and functions can be recognized more easily.

Attention

Argument-Dependent Name Lookup (Koenig Lookup)

Using the namespace prefix seqan:: is not really necessary in all places.
In many cases, the Koenig lookup rule in C++ for functions makes this unnecessary.
Consider the following, compiling, example.

seqan::String<char>s="example";unsignedi=length(s);

Here, the function length does not have a namespace prefix.
The code compiles nevertheless.
The compiler automatically looks for a function length in the namespace of its arguments.

Note that we follow the rules for variable, function, and class names as outlined in the SeqAn style guide.
For example:
1. variables and functions use lower case,
2. struct, enum and classes use CamelCase,
3. metafunctions start with a capital letter, and
4. metafunction values are UPPERCASE.

Depending on your operating system you have different alternatives to create a demo application.
An in depth description can be found in GettingStarted.

Solution

Click ‚Äė‚Äômore‚Ä¶‚Äô‚Äė

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>usingnamespaceseqan;intmain(){// InitializationString<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score;resize(score,length(text)-length(pattern)+1);// Computation of the similarities// Iteration over the text (outer loop)for(unsignedi=0;i<length(text)-length(pattern)+1;++i){intlocalScore=0;// Iteration over the pattern for character comparisonfor(unsignedj=0;j<length(pattern);++j){if(text[i+j]==pattern[j])++localScore;}score[i]=localScore;}// Printing the resultfor(unsignedi=0;i<length(score);++i)std::cout<<score[i]<<" ";std::cout<<std::endl;return0;}

The String class is one of the most fundamental classes in SeqAn, which comes as no surprise since SeqAn is used to analyse sequences (there is an extra tutorial for SeqAn sequences and alphabets).

In contrast to the popular string classes of Java or C++, SeqAn provides different string implementations and different alphabets for its strings.
There is one string implementation that stores characters in memory, just like normal C++ strings.
Another string implementation stores the characters on disk and only keeps a part of the sequence in memory.
For alphabets, you can use strings of nucleotides, such as genomes, or you can use strings of amino acids, for example.

SeqAn uses template functions and template classes to implement the different types of strings using the generic programming paradigm.
Template functions/classes are normal functions/classes with the additional feature that one passes the type of a variable as well as its value (see also: templates in cpp).
This means that SeqAn algorithms and data structures are implemented in such a way that they work on all types implementing an informal interface (see information box below for more details).
This is similar to the philosophy employed in the C++ STL (Standard Template Library).

The following two lines make use of template programming to define two strings of type char, a text and a pattern.

// InitializationString<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";

In order to store the similarities between the pattern and different text positions we additionally create a string storing integer values.

String<int>score;

Note that in contrast to the first two string definitions we do not know the values of the different positions in the string in advance.
In order to dynamically adjust the length of the new string to the text we can use the function resize.
The resize function is not a member function of the string class because SeqAn is not object oriented in the typical sence (we will see later how we adapt SeqAn to object oriented programming).
Therefore, instead of writing string.resize(newLength) we use resize(string,newLength).

resize(score,length(text)-length(pattern)+1);

Note

Global function interfaces.

SeqAn uses global interfaces for its data types/classes.
Generally, you have to use function(variable) instead of variable.function().

This has the advantage that we can extend the interface of a type outside of its definition.
For example, we can provide a length() function for STL containers std::string<T> and std::vector<T> outside their class files.
We can use such global functions to make one data type have the same interface as a second.
This is called adaption.

Additionally, we can use one function definition for several data types.
For example, the alignment algorithms in SeqAn are written such that we can compute alignments using any String with any alphabet:
There are more than 5 String variants in SeqAn and more than 8 built-in alphabets.
Thus, one implementation can be used for more than 40 different data types!

After the string initializations it is now time for the similarity computation.
In this toy example we simply take the pattern and shift it over the text from left to right.
After each step, we check how many characters are equal between the corresponding substring of the text and the pattern.
We implement this using two loops; the outer one iterates over the given text and the inner loop over the given pattern:

// Computation of the similarities// Iteration over the text (outer loop)for(unsignedi=0;i<length(text)-length(pattern)+1;++i){intlocalScore=0;// Iteration over the pattern for character comparisonfor(unsignedj=0;j<length(pattern);++j){if(text[i+j]==pattern[j])++localScore;}score[i]=localScore;}

There are two things worth mentioning here: (1) SeqAn containers or strings start at position 0 and (2) you will notice that we use ++variable instead of variable++ wherever possible.
The reason is that ++variable is slightly faster than its alternative, since the alternative needs to make a copy of itself before returning the result.

In the last step we simply print the result that we stored in the variable ````score on screen.
This gives the similarity of the pattern to the string at each position.

At this point, we have already created a working solution!
However, in order to make it easier to maintain and reuse parts of the code we need to export them into functions.
In this example the interesting piece of code is the similarity computation, which consists of an outer and inner loop.
We encapsulate the outer loop in function computeScore and the inner loop in function computeLocalScore as can be seen in the following code.

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>usingnamespaceseqan;intcomputeLocalScore(String<char>subText,String<char>pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}String<int>computeScore(String<char>text,String<char>pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);for(unsignedi=0;i<length(score);++i)std::cout<<score[i]<<" ";std::cout<<std::endl;return0;}

The function computeScore() now contains the fundamental part of the code and can be reused by other functions.
The input arguments are two strings.
One is the pattern itself and one is a substring of the text.
In order to obtain the substring we can use the function infix implemented in SeqAn.
The function call infix(text,i,j) generates a substring equal to text[i...j-1], e.g. infix(text,1,5) equals ‚Äúello‚ÄĚ, where text is ‚ÄúHello World‚ÄĚ.
To be more precise, infix() generates a Infix which can be used as a string, but is implemented using pointers such that no copying is necessary and running time and memory is saved.

Replace the code in your current file by the code above and encapsulate the print instructions.

Hint

The function head should look like this:

voidprint(String<int>text)

Solution

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>usingnamespaceseqan;intcomputeLocalScore(String<char>subText,String<char>pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}String<int>computeScore(String<char>text,String<char>pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}voidprint(String<int>text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(score);return0;}

Both the text and the pattern are passed by value.
This means that both the text and the pattern are copied when the function is called, which consumes twice the memory.
This can become a real bottleneck since copying longer sequences is very memory and time consuming, think of the human genome, for example.

Instead of copying we could use references.
A reference in C++ is created using an ampersand sign (&) and creates an alias to the referenced value.
Basically, a reference is a pointer to an object which can be used just like the referenced object itself.
This means that when you change something in the reference you also change the original object it came from.
But there is a solution to circumvent this modification problem as well, namely the word const.
A const object cannot be modified.

Important

If an object does not need to be modified make it an nonmodifiably object using the keyword const.
This makes it impossible to unwillingly change objects, which can be really hard to debug.
Therefore it is recommended to use it as often as possible.

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>usingnamespaceseqan;intcomputeLocalScore(String<char>const&subText,String<char>const&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}String<int>computeScore(String<char>const&text,String<char>const&pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}voidprint(String<int>const&text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(score);return0;}

As mentioned earlier, there is another issue: the function computeScore only works for Strings having the alphabet char.
If we wanted to use it for Dna or AminoAcid strings then we would have to reimplement it even though the only difference is the signature of the function.
All used functions inside computeScore can already handle the other datatypes.

The more appropriate solution is a generic design using templates, as often used in the SeqAn library.
Instead of specifying the input arguments to be references of strings of char s we could use references of template arguments as shown in the following lines:

The first line above specifies that we create a template function with two template arguments TText and TPattern.
At compile time the template arguments are then replace with the correct types.
If this line was missing the compiler would expect that there are types TText and TPattern with definitions.

Now the function signature is better in terms of memory consumption, time efficiency, and generality.

Important

The SeqAn Style Guide

The SeqAn style guide gives rules for formatting and structuring C++ code as well as naming conventions.
Such rules make the code more consistent, easier to read, and also easier to use.

Naming Scheme.
Variable and function names are written in lowerCamelCase, type names are written in UpperCamelCase.
Constants and enum values are written in UPPER_CASE.
Template variable names always start with ‚ÄėT‚Äô.

Function Parameter Order.
The order is (1) output, (2) non-const input (e.g. file handles), (3) input, (4) tags.
Output and non-const input can be modified, the rest is left untouched and either passed by copy or by const-reference (const&).

Global Functions.
With the exception of constructors and a few operators that have to be defined in-class, the interfaces in SeqAn use global functions.

No Exceptions.
The SeqAn interfaces do not throw any exceptions.

While we are trying to make the interfaces consistent with our style guide, some functions have incorrect parameter order.
This will change in the near future to be more in line with the style guide.

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>#include<seqan/score.h>usingnamespaceseqan;template<typenameTText,typenameTPattern>intcomputeLocalScore(TTextconst&subText,TPatternconst&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}template<typenameTText,typenameTPattern>String<int>computeScore(TTextconst&text,TPatternconst&pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}voidprint(String<int>const&text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(score);return0;}

There is another huge advantage of using templates: we can specialize a function without touching the existing function.
In our working example it might be more appropriate to treat AminoAcid sequences differently.
As you probably know, there is a similarity relation on amino acids: Certain amino acids are more similar to each other, than others.
Therefore we want to score different kinds of mismatches differently.
In order to take this into consideration we simple write a computeLocalScore() function for AminoAcid strings.
In the future whenever ‚ÄėcomputerScore‚Äô is called always the version above is used unless the second argument is of type String-AminoAcid.
Note that the second template argument was removed since we are using the specific type String-AminoAcid.

In order to score a mismatch we use the function score() from the SeqAn library.
Note that we use the Blosum62 matrix as a similarity measure.
When looking into the documentation of score you will notice that the score function requires a argument of type Score.
This object tells the function how to compare two letters and there are several types of scoring schemes available in SeqAn (of course, you can extend this with your own).
In addition, because they are so frequently used there are shortcuts as well.
For example Blosum62 is really a shortcut for Score<int,ScoreMatrix<AminoAcid,Blosum62_>>, which is obviously very helpful.
Other shortcuts are DnaString for String<Dna> (sequence tutorial), CharString for String<char>, ‚Ä¶

Tip

Template Subclassing

The main idea of template subclassing is to exploit the C++ template matching mechanism.
For example, in the following code, the function calls (1) and (3) will call the function myFunction() in variant (A) while the function call (2) will call variant (B).

Provide a generic print function which is used when the input type is not String<int>.

Hint

Keep your current implementation and add a second function.
Don‚Äôt forget to make both template functions.
Include <seqan/score.h> as well.

Solution

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>#include<seqan/score.h>usingnamespaceseqan;template<typenameTText>intcomputeLocalScore(TTextconst&subText,String<AminoAcid>const&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)localScore+=score(Blosum62(),subText[i],pattern[i]);returnlocalScore;}template<typenameTText,typenameTPattern>intcomputeLocalScore(TTextconst&subText,TPatternconst&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}template<typenameTText,typenameTPattern>String<int>computeScore(TTextconst&text,TPatternconst&pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}template<typenameTText>voidprint(TTextconst&text){std::cout<<text<<std::endl;}voidprint(String<int>const&text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(text);print(pattern);print(score);return0;}

Having a closer look you will notice that there is a default constructor call (MyersHirschberg() ) within a function call.
Using this mechanism one can specify which function to call at compile time.
The MyersHirschberg() `` is only a tag to determine which specialisation of the globalAligment function to call.

If you want more information on tags then read on otherwise you are now ready to explore SeqAn in more detail and continue with one of the other tutorials.

There is another use case of templates and function specialization.

This might be useful in a print() function, for example.
In some scenarios, we only want to print the position where the maximal similarity between pattern and text is found.
In other cases, we might want to print the similarities of all positions.
In SeqAn, we use tag-based dispatching to realize this.
Here, the type of the tag holds the specialization information.

Tip

Tag-Based Dispatching

You will often see tags in SeqAn code, e.g. Standard().
These are parameters to functions that are passed as const-references.
They are not passed for their values but for their type only.
This way, we can select different specializations at compile time in a way that plays nicely together with metafunctions, template specializations, and an advanced technique called [[Tutorial/BasicTechniques| metaprogramming]].

If we call print() with something different than MaxOnly then we print all the positions with their similarity, because the generic template function accepts anything as the template argument.
On the other hand, if we call print with MaxOnly only the positions with the maximum similarity as well as the maximal similarity will be shown.

// Provide a print function that prints pairs of positions and their score if the score is greater than 0.#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>#include<seqan/score.h>usingnamespaceseqan;template<typenameTText>intcomputeLocalScore(TTextconst&subText,String<AminoAcid>const&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)localScore+=score(Blosum62(),subText[i],pattern[i]);returnlocalScore;}template<typenameTText,typenameTPattern>intcomputeLocalScore(TTextconst&subText,TPatternconst&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}template<typenameTText,typenameTPattern>String<int>computeScore(TTextconst&text,TPatternconst&pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}template<typenameTText>voidprint(TTextconst&text){std::cout<<text<<std::endl;}voidprint(String<int>const&text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}template<typenameTText,typenameTSpec>voidprint(TTextconst&text,TSpecconst&/*tag*/){print(text);}structMaxOnly{};template<typenameTText>voidprint(TTextconst&score,MaxOnlyconst&/*tag*/){intmaxScore=score[0];String<int>output;appendValue(output,0);for(unsignedi=1;i<length(score);++i){if(score[i]>maxScore){maxScore=score[i];clear(output);resize(output,1,i);}elseif(score[i]==maxScore)appendValue(output,i);}print(output);}structGreaterZero{};template<typenameTText>voidprint(TTextconst&score,GreaterZeroconst&/*tag*/){String<Pair<int>>output;for(unsignedi=1;i<length(score);++i)if(score[i]>0)appendValue(output,Pair<int>(i,score[i]));for(unsignedi=0;i<length(output);++i)std::cout<<"("<<output[i].i1<<"; "<<output[i].i2<<") ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(text);print(pattern);print(score);print(score,MaxOnly());print(score,GreaterZero());// And now for a protein patternString<AminoAcid>protein="tutorial";String<int>proteinScore=computeScore(text,protein);print(text);print(protein);print(proteinScore);print(proteinScore,MaxOnly());print(proteinScore,GreaterZero());return0;}

Obviously this is only a toy example in which we could have named the two print() functions differently.
However, often this is not the case when the programs become more complex.
Because SeqAn is very generic we do not know the datatypes of template functions in advance.
This would pose a problem because the function call of function b() in function a() may depend on the data types of the template arguments of function a().

Don‚Äôt worry if you have not fully understood the last section.
If you have ‚Äď perfect.
In any case the take home message is that you use data types for class specializations and if you see a line of code in which the default constructor is written in a function call this typical means that the data type is important to distinct between different function implementations.

Now you are ready to explore more of the SeqAn library.
There are several tutorials which will teach you how to use the different SeqAn data structures and algorithms.
Below you find the complete code for our example with the corresponding output.

#include<iostream>#include<seqan/file.h>#include<seqan/sequence.h>#include<seqan/score.h>usingnamespaceseqan;template<typenameTText>intcomputeLocalScore(TTextconst&subText,String<AminoAcid>const&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)localScore+=score(Blosum62(),subText[i],pattern[i]);returnlocalScore;}template<typenameTText,typenameTPattern>intcomputeLocalScore(TTextconst&subText,TPatternconst&pattern){intlocalScore=0;for(unsignedi=0;i<length(pattern);++i)if(subText[i]==pattern[i])++localScore;returnlocalScore;}template<typenameTText,typenameTPattern>String<int>computeScore(TTextconst&text,TPatternconst&pattern){String<int>score;resize(score,length(text)-length(pattern)+1,0);for(unsignedi=0;i<length(text)-length(pattern)+1;++i)score[i]=computeLocalScore(infix(text,i,i+length(pattern)),pattern);returnscore;}template<typenameTText>voidprint(TTextconst&text){std::cout<<text<<std::endl;}voidprint(String<int>const&text){for(unsignedi=0;i<length(text);++i)std::cout<<text[i]<<" ";std::cout<<std::endl;}template<typenameTText,typenameTSpec>voidprint(TTextconst&text,TSpecconst&/*tag*/){print(text);}structMaxOnly{};template<typenameTText>voidprint(TTextconst&score,MaxOnlyconst&/*tag*/){intmaxScore=score[0];String<int>output;appendValue(output,0);for(unsignedi=1;i<length(score);++i){if(score[i]>maxScore){maxScore=score[i];clear(output);resize(output,1,i);}elseif(score[i]==maxScore)appendValue(output,i);}print(output);}structGreaterZero{};template<typenameTText>voidprint(TTextconst&score,GreaterZeroconst&/*tag*/){String<Pair<int>>output;for(unsignedi=1;i<length(score);++i)if(score[i]>0)appendValue(output,Pair<int>(i,score[i]));for(unsignedi=0;i<length(output);++i)std::cout<<"("<<output[i].i1<<"; "<<output[i].i2<<") ";std::cout<<std::endl;}intmain(){String<char>text="This is an awesome tutorial to get to know SeqAn!";String<char>pattern="tutorial";String<int>score=computeScore(text,pattern);print(text);print(pattern);print(score);print(score,MaxOnly());print(score,GreaterZero());// And now for a protein patternString<AminoAcid>protein="tutorial";String<int>proteinScore=computeScore(text,protein);print(text);print(protein);print(proteinScore);print(proteinScore,MaxOnly());print(proteinScore,GreaterZero());return0;}

You will learn how to use the ArgumentParser class to parse command line arguments.
This tutorial is a walk-through with links into the API documentation and also meant as a source for copy-and-paste code.

The following small program will (1) setup a ArgumentParser object named parser, (2) parse the command line, (3) exit the program if there were errors or the user requested a functionality that is already built into the command line parser, and (4) printing the settings given from the command line.
Such functionality is printing the help, for example.

#include<iostream>#include<seqan/arg_parse.h>intmain(intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;// Extract option values and print them.unsignedperiod=0;getOptionValue(period,parser,"period");booltoUppercase=isSet(parser,"uppercase");seqan::CharStringtext;getArgumentValue(text,parser,0);std::cout<<"period \t"<<period<<'\n'<<"uppercase\t"<<toUppercase<<'\n'<<"text \t"<<text<<'\n';return0;}

Let us first play a bit around with the program before looking at it in detail.

Let us look at this program in detail now. The required SeqAn module is seqan/arg_parse.h.
After inclusion, we can create an ArgumentParser object:

// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");

Then, we define a positional argument using the function addArgument.
The function accepts the parser and an ArgParseArgument object.
We call the ArgParseArgument constructor with two parameters: the type of the argument (a string), and a label for the documentation.

addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));

The ArgParseOption constructor is called in two different variants.
Within the first addOption call, we construct an integer option with a short and long name, a documentation string, and give it the label ‚ÄúINT‚ÄĚ.
The second option is a flag (indicated by not giving a type) with a short and a long name and a description.

We then check the result of the parsing operation.
The result is seqan::ArgumentParser::PARSE_ERROR if there was a problem with the parsing.
Otherwise, it is seqan::ArgumentParser::PARSE_OK if there was no problem and no special functionality of the argument parser was triggered.
The command line parser automatically adds some arguments, such as --help.
If such built-in functionality is triggered, it will return a value that is neither PARSE_ERROR nor PARSE_OK.

The following two lines have the following behaviour.
If the parsing went through and no special functionality was triggered then the branch is not taken.
Otherwise, the method main() is left with 1 in case of errors and with 0 in case special behaviour was triggered (e.g. the help was printed).

Finally, we access the values from the command line using the ArgumentParser.
The function getOptionValue allows us to access the values from the command line after casting into C++ types.
The function isSet allows us to query whether a given argument was set on the command line.

You have to mark an option to be a list if you want to be able to collect multiple values for it from the command line.
Consider the following program call:

# program -a 1 -a 2 -a 3

If the option a is not a list then the occurence -a3 overwrites all previous settings.

However, if a is marked to be a list, then all values (1, 2, and 3) are stored as its values.
We can get the number of elements using the function getOptionValueCount and then access the individual arguments using the function getOptionValue.
You can mark an option and arguments to be lists by using the isList parameter to the ArgParseArgument and ArgParseOption constructors.

For arguments, only the first or the last argument or none can be a list but not both.
Consider this program call:

# program arg0 arg1 arg2 arg3

For example, if the program has three arguments and the first one is a list then arg0 and arg1 would be the content of the first argument.
If it has two arguments and the last one is a list then arg1, arg2, and arg3 would be the content of the last argument.

Adjust the program from above to also accept an option to convert characters to lower case, just as it accepts options to convert characters to upper case.
The long name should be --lowercase, the short name should be -L.
As for the --uppercase option, the program should print whether the flag was set or not.

Hint

Copy the two lines for defining the --uppercase option and replace the strings appropriately.

Solution

#include<iostream>#include<seqan/arg_parse.h>intmain(intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;// Extract option values and print them.unsignedperiod=0;getOptionValue(period,parser,"period");booltoUppercase=isSet(parser,"uppercase");booltoLowercase=isSet(parser,"lowercase");seqan::CharStringtext;getArgumentValue(text,parser,0);std::cout<<"period \t"<<period<<'\n'<<"uppercase\t"<<toUppercase<<'\n'<<"lowercase\t"<<toLowercase<<'\n'<<"text \t"<<text<<'\n';return0;}

Adjust the previous program to accept default values by adding the setDefaultValue() line from above into your program.

Solution

#include<iostream>#include<seqan/arg_parse.h>intmain(intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;// Extract option values and print them.unsignedperiod=0;getOptionValue(period,parser,"period");booltoUppercase=isSet(parser,"uppercase");booltoLowercase=isSet(parser,"lowercase");seqan::CharStringtext;getArgumentValue(text,parser,0);std::cout<<"period \t"<<period<<'\n'<<"uppercase\t"<<toUppercase<<'\n'<<"lowercase\t"<<toLowercase<<'\n'<<"text \t"<<text<<'\n';return0;}

Instead of just printing the options back to the user, we should actually store them.
To follow best practice, we should not use global variables for this but instead pass them as parameters.

We will thus create a ModifyStringOptions struct that encapsulates the settings the user can give to the modify_string program.
Note that we initialize the variables of the struct with initializer lists, as it is best practice in modern C++.

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),toUppercase(false),toLowercase(false){}};intmain(intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;// Extract option values and print them.ModifyStringOptionsoptions;getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");getArgumentValue(options.text,parser,0);std::cout<<"period \t"<<options.period<<'\n'<<"uppercase\t"<<options.toUppercase<<'\n'<<"lowercase\t"<<options.toLowercase<<'\n'<<"text \t"<<options.text<<'\n';return0;}

As a next step towards a cleaner program, we should extract the argument parsing into its own function, e.g. call it parseCommandLine().
Following the style guide (C++ Code Style), we first pass the output parameter, then the input parameters.
The return value of our function is a seqan::ArgumentParser::ParseResult such that we can differentiate whether the program can go on, the help was printed and the program is to exit with success, or there was a problem with the passed argument and the program is to exit with an error code.

Also, note that we should check that the user cannot specify both to-lowercase and to-uppercase.
This check cannot be performed by the ArgumentParser by itself but we can easily add this check.
We add this functionality to the parseCommandLine() function.

Click more‚Ä¶ to see the updated program.

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// We require one argument.addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));// Define OptionsaddOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values. getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");getArgumentValue(options.text,parser,0);// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::cout<<"period \t"<<options.period<<'\n'<<"uppercase\t"<<options.toUppercase<<'\n'<<"lowercase\t"<<options.toLowercase<<'\n'<<"text \t"<<options.text<<'\n';return0;}

The command line parsing part of our program is done now.
Let us now add a function modifyString() that is given a ModifyStringOptions object and text and modifies the text.
We simply use the C standard library functios toupper() and tolower() from the header <cctype> for converting to upper and lower case.

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// We require one argument.addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));// Define OptionsaddOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values. getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");getArgumentValue(options.text,parser,0);// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}seqan::CharStringmodifyString(seqan::CharStringconst&text,ModifyStringOptionsconst&options){seqan::CharStringresult;if(options.toLowercase){for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,tolower(text[i]));elseappendValue(result,text[i]);}}else{for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,toupper(text[i]));elseappendValue(result,text[i]);}}returnresult;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::cout<<modifyString(options.text,options)<<'\n';return0;}

Use the function setMinValue to set a minimal value of 1 for the parameter --period.

Solution

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;unsignedrangeBegin,rangeEnd;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),rangeBegin(0),rangeEnd(0),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// We require one argument.addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));// Define OptionsaddOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setMinValue(parser,"period","1");setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("r","range","Range of the text to modify.",seqan::ArgParseArgument::INTEGER,"INT",false,2));addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values.getOptionValue(options.period,parser,"period");getOptionValue(options.rangeBegin,parser,"range",0);getOptionValue(options.rangeEnd,parser,"range",1);options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");seqan::getArgumentValue(options.text,parser,0);// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}seqan::CharStringmodifyString(seqan::CharStringconst&text,ModifyStringOptionsconst&options){seqan::CharStringresult;if(options.toLowercase){for(unsignedi=0;i<length(text);++i){if(i>=options.rangeBegin&&i<options.rangeEnd&&(i%options.period==0u))appendValue(result,tolower(text[i]));elseappendValue(result,text[i]);}}else{for(unsignedi=0;i<length(text);++i){if(i>=options.rangeBegin&&i<options.rangeEnd&&(i%options.period==0u))appendValue(result,toupper(text[i]));elseappendValue(result,text[i]);}}returnresult;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::cout<<modifyString(options.text,options)<<'\n';return0;}

Sometimes, it is useful to give a list of valid values for a command line option.
You can give it as a space-separated list in a string to setValidValues.
The check whether the value from the command line is valid is case sensitive.

We could use ArgParseArgument::STRING to specify input and output files.
However, there are two special argument/option types ArgParseArgument::INPUT_FILE and ArgParseArgument::OUTPUT_FILE that are more suitable:

In the near future, we plan to add basic checks for whether input files exist and are readable by the user.
You will still have to check whether opening was successful when actually doing this but the program will fail earlier if the source file or target location are not accessible.
The user will not have to wait for the program to run through to see that he mistyped the output directory name, for example, and you do not have to write this check.

For workflow engine integration, the input and output file options and arguments will be converted into appropriate input and output ports of the nodes.

You can use the previously introduced restrictions to specify what kind of files you expect and the ArgumentParser will check while parsing if the correct file type was provided.

Here is an example for defining input and output file arguments:

seqan::addOption(parser,seqan::ArgParseOption("I","input-file","Path to the input file",seqan::ArgParseArgument::INPUT_FILE,"IN"));seqan::addOption(parser,seqan::ArgParseOption("O","output-file","Path to the output file",seqan::ArgParseArgument::OUTPUT_FILE,"OUT"));

Replace the argument TEXT by a command line option -I/--input-file in the program above.
The program should then read in the text instead of using the command line argument.

Hint

We will also replace the text member of ModifyStringOptions, you might wish to do the same.

Solution

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;unsignedrangeBegin,rangeEnd;booltoUppercase;booltoLowercase;seqan::CharStringinputFileName;ModifyStringOptions():period(1),rangeBegin(0),rangeEnd(0),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// Define OptionsaddOption(parser,seqan::ArgParseOption("I","input-file","A text file that will printed with the modifications applied.",seqan::ArgParseArgument::INPUT_FILE));setValidValues(parser,"input-file","txt");setRequired(parser,"input-file");addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setMinValue(parser,"period","1");setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values.getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");getOptionValue(options.inputFileName,parser,"input-file");// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}seqan::CharStringmodifyString(seqan::CharStringconst&text,ModifyStringOptionsconst&options){seqan::CharStringresult;if(options.toLowercase){for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,tolower(text[i]));elseappendValue(result,text[i]);}}else{for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,toupper(text[i]));elseappendValue(result,text[i]);}}returnresult;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::fstreaminFile(toCString(options.inputFileName),std::ios::binary|std::ios::in);if(inFile.good()){std::cerr<<"ERROR: Could not open input file "<<options.inputFileName<<'\n';return1;}seqan::CharStringtext;while(inFile.good()){charc=inFile.get();if(inFile.good())appendValue(text,c);}std::cout<<modifyString(text,options);return0;}

Add a command line option --range to the ArgumentParser in the program above.
Modify the function modifyString() such that only parameters in the given range are changed.

Hint

We will add two unsigned members rangeBegin and rangeEnd to the ModifyStringOptions struct, you might wish to do the same.

Solution

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;unsignedrangeBegin,rangeEnd;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),rangeBegin(0),rangeEnd(0),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// We require one argument.addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));// Define OptionsaddOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setMinValue(parser,"period","1");setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values.getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");seqan::getArgumentValue(options.text,parser,0);// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}seqan::CharStringmodifyString(seqan::CharStringconst&text,ModifyStringOptionsconst&options){seqan::CharStringresult;if(options.toLowercase){for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,tolower(text[i]));elseappendValue(result,text[i]);}}else{for(unsignedi=0;i<length(text);++i){if(i%options.period==0u)appendValue(result,toupper(text[i]));elseappendValue(result,text[i]);}}returnresult;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::cout<<modifyString(options.text,options)<<'\n';return0;}

Another very useful feature of ArgumentParser is that you can embed rich documentation into your programs.
You can set the short description, the version string, date, synopsis and add text documentation settings.

Let us first set the short description, version string, and date in our program from above.
We insert the following lines just after the declaration of the variable parser.

The formatting of command line parameters might seem strange, at first:
Font operators start with \f (which means that they start with "\\f" in in C++ string literals).
The \\f is followed by the format specifier.
The format specifier can be one of I, B, and P.
I selects italic text (underlined on the shell), B selects bold and P resets the formatting to normal text.
These font operators are legacies of man pages from Unix and offered a simple-to-implement solution to text formatting.

For example, "Words\\fBwere\\fPmadefor\\fIbeing\\fPwritten!" would result in the formatted string ‚ÄúWords were made for being written!‚ÄĚ.

Note that formatting the command line relies on ANSI escape codes which is not supported by modern Windows versions.
If you are using Windows, you will not see bold or underlined text.

The argument parser will add some options of its own, for example for printing the help and displaying version information.
To separate our arguments from the autogenerated ones, we add the following line.
This line will introduce the section ‚ÄúModification Options‚ÄĚ in the Description section of the output.

seqan::addSection(parser,"Modification Options");

Finally, we will add a section with examples.
Add the following lines just before the line with the parse() function call.

seqan::addTextSection(parser,"Examples");seqan::addListItem(parser,"\\fBmodify_string\\fP \\fB-U\\fP \\fIveryverylongword\\fP","Print upper case version of \"veryverylongword\"");seqan::addListItem(parser,"\\fBmodify_string\\fP \\fB-L\\fP \\fB-i\\fP \\fI3\\fP \\fIveryverylongword\\fP","Print \"veryverylongword\" with every third character ""converted to upper case.");

That were a lot of changes!
Click more‚Ä¶ to see the complete program.

#include<iostream>#include<seqan/arg_parse.h>structModifyStringOptions{unsignedperiod;booltoUppercase;booltoLowercase;seqan::CharStringtext;ModifyStringOptions():period(1),toUppercase(false),toLowercase(false){}};seqan::ArgumentParser::ParseResultparseCommandLine(ModifyStringOptions&options,intargc,charconst**argv){// Setup ArgumentParser.seqan::ArgumentParserparser("modify_string");// Set short description, version, and date.setShortDescription(parser,"String Modifier");setVersion(parser,"1.0");setDate(parser,"July 2012");// Define usage line and long description.addUsageLine(parser,"[\\fIOPTIONS\\fP] \"\\fITEXT\\fP\"");addDescription(parser,"This program allows simple character modifications to ""each i-th character.");// We require one argument.addArgument(parser,seqan::ArgParseArgument(seqan::ArgParseArgument::STRING,"TEXT"));// Define Options -- Section Modification OptionsaddSection(parser,"Modification Options");addOption(parser,seqan::ArgParseOption("i","period","Period to use for the index.",seqan::ArgParseArgument::INTEGER,"INT"));setDefaultValue(parser,"period","1");addOption(parser,seqan::ArgParseOption("U","uppercase","Select to-uppercase as operation."));addOption(parser,seqan::ArgParseOption("L","lowercase","Select to-lowercase as operation."));// Add Examples Section.addTextSection(parser,"Examples");addListItem(parser,"\\fBmodify_string\\fP \\fB-U\\fP \\fIveryverylongword\\fP","Print upper case version of \"veryverylongword\"");addListItem(parser,"\\fBmodify_string\\fP \\fB-L\\fP \\fB-i\\fP \\fI3\\fP ""\\fIveryverylongword\\fP","Print \"veryverylongword\" with every third character ""converted to upper case.");// Parse command line.seqan::ArgumentParser::ParseResultres=seqan::parse(parser,argc,argv);// Only extract options if the program will continue after parseCommandLine()if(res!=seqan::ArgumentParser::PARSE_OK)returnres;// Extract option values.getOptionValue(options.period,parser,"period");options.toUppercase=isSet(parser,"uppercase");options.toLowercase=isSet(parser,"lowercase");seqan::getArgumentValue(options.text,parser,0);// If both to-uppercase and to-lowercase were selected then this is an error.if(options.toUppercase&&options.toLowercase){std::cerr<<"ERROR: You cannot specify both to-uppercase and to-lowercase!\n";returnseqan::ArgumentParser::PARSE_ERROR;}returnseqan::ArgumentParser::PARSE_OK;}seqan::CharStringmodifyString(seqan::CharStringconst&text,ModifyStringOptionsconst&options){seqan::CharStringresult;if(options.toLowercase){for(unsignedi=0;i<length(text);++i)appendValue(result,tolower(text[i]));}else{for(unsignedi=0;i<length(text);++i)appendValue(result,toupper(text[i]));}returnresult;}intmain(intargc,charconst**argv){// Parse the command line.ModifyStringOptionsoptions;seqan::ArgumentParser::ParseResultres=parseCommandLine(options,argc,argv);// If parsing was not successful then exit with code 1 if there were errors.// Otherwise, exit with code 0 (e.g. help was printed).if(res!=seqan::ArgumentParser::PARSE_OK)returnres==seqan::ArgumentParser::PARSE_ERROR;std::cout<<modifyString(options.text,options)<<'\n';return0;}

Let us look at the resulting documentation.
Simply call the new program with the --help option.

Also, there is an undocumented option called --export-help that is automatically added by ArgumentParser.
You can call it with the values html and man.
If the option is set then the argument parser will print the documentation as HTML or man format (man pages are a widely used format for Unix documentation).

HTML can be displayed by any web browser, man pages can be displayed using the program man.
Note that when opening a file using man, you have to give the file name either as an absolute or a relative path.
Otherwise, it would try to look up the topic modify_string.man.
To view the generated man page use:

# man ./modify_string.man

Below, you can see a part of the rendered HTML and man pages generated by the commands above.

With the seqan-2.3.0 release applications, using the ArgumentParser, check SeqAn servers for version updates. The functionality helps getting new versions out to users faster. It is also used to inform application developers of new versions of the SeqAn library which means that applications ship with less bugs.

The version information you receive depends on whether you are an application user or developer.
We differentiate this by inquiring the NDEBUG (no debug) macro.

#. Case: NDEBUGis set. Its the default in our application and represents that you ara a user.
The only message you will eventually encounter is the following:

[APP INFO] :: There is a newer version of this application available.[APP INFO] :: If this app is developed by SeqAn, visit www.seqan.de for updates.[APP INFO] :: If you don't want to recieve this message again set --version_check OFF

#. Case: NDEBUGis NOT set. If you build one of our application or your own one in debug mode, we will consider you as a developer.
Therefore we will inform you whenever a new library version is available:

[SEQAN INFO] :: There is a newer SeqAn version available! [SEQAN INFO] :: Please visit www.seqan.de for an update or inform the developer of this app. [SEQAN INFO] :: If you don't want to recieve this message again set --version-check OFFIf you are working on your own application, using the SeqAn ArgumentParser, we will inform you about the possibility to register your application with us. This will make the distribution of your application version simple and convenient.

[SEQAN INFO] :: Thank you for using SeqAn![SEQAN INFO] :: You might want to regsiter you app for support and version check features?[SEQAN INFO] :: Just send us an email to seqan@team.fu-berlin.de with your app name and version number.[SEQAN INFO] :: If you don't want to recieve this message anymore set --version_check OFF

The process of checking for a new version happens at most once a day and takes at most three seconds enforced by an internal timeout.

Note

The runtime of your application might be slightly affected by the process of checkng the version.
You might want to temporarily switch off the option while doing sensible performance measurements (--version-checkOFF).

The following information is transmitted to the servers solely via the URL:

The application name and its version

The SeqAn version that the application was built with

The operating system type (Linux, macOS, Windows or BSD)

The CPU type (32 or 64bit)

The purpose of this transmission is to provide accurate update data for your app.
Beyond the update information, we may count the total number of version requests and may also resolve them to geographical regions.
This may be used for anonymized analysis by the SeqAn team, but raw data is never shared with third parties.

Attention

There is no form of user identification and no tracking.
IP-Addresses are never stored permanently.
SeqAn collects no information regarding your use of the application, like selected options or arguments provided, and of course no information on your files!

You will learn about the SeqAn sequence concept and its main class String as well as the class Segment.
After completing this tutorial, you will be able to use important functionalities of sequences in SeqAn and you will be ready to continue with the more specific tutorials, e.g. Alignment, or Pairwise Sequence Alignment.

Sequences are the core concept of SeqAn.
A sequence is a container that stores an ordered list of values.
In SeqAn, there are three kinds of sequences: Strings, Sequence Adaptions and Segments.

The String class is one of the most fundamental classes in SeqAn.
It is designed as a generic data structure that can be instantiated for all kinds of values, both simple (e.g. char, Dna, AminoAcid) and non-simple value types (e.g. Tuple, String).
With sequence adaptions, SeqAn offers an interface for accessing data types that are not part of SeqAn, namely standard library strings and c-style char arrays.
Thus those built-in types can be handled in a similar way as SeqAn strings, for example with the length function.
Segments are contiguous subsequences that represent parts of other sequences.

This tutorial will deal with the SeqAn sequence classes String and Segment.

Let‚Äôs first have a look at a simple example on how to define a String.
The type of the contained value is specified by the first template argument, e.g. char or int.

String<char>myText;// A string of characters.String<int>myNumbers;// A string of integers.

To fill the string with contents, we can simply assign a string literal to the created variable:

myText="Hello SeqAn!";

Any type that provides a default constructor, a copy constructor and an assignment operator can be used as the alphabet / contained type of a String.
This includes the C++ POD types, e.g. char, int, double etc., or even more complex types complex types, such as Strings.

String<String<char>>myStringList;// A string of character strings.

Hint

Nested Sequences (aka ‚ÄúStrings of Strings‚ÄĚ)

A collection of sequences can either be stored in a sequence of sequences, for example in a String<String<char>>, or in a StringSet.
The latter one allows for more auxiliary functionalities to improve the efficiency of working with large sequence collections.
You can learn more about it in the tutorial String Sets.

The SeqAn String implementation provides the common C++ operators that you know already from the vector class of the STL.
For example:

String<Dna>dnaSeq="TATA";dnaSeq+="CGCG";std::cout<<dnaSeq<<std::endl;

TATACGCG

Each sequence object has a capacity, i.e. the maximum length of a sequence that can be stored in this object.
While some sequence types have a fixed capacity, the capacity of other sequence classes like Alloc String or std::basic_string can be changed at runtime.
The capacity can be set explicitly by functions such as reserve or resize.
It can also be set implicitly by functions like append or replace, if the operation‚Äôs result exceeds the length of the target string.

In the following example we create a String of Dna5String. We first set the new length of the container with resize to two elements.
After assigning two elements we append one more element with appendValue.
In the last step the capacity is implicitly changed.

// And to check if your output is correct,// use the given SeqAn function reverseComplement(),// which modifies the sequence in-place:reverseComplement(genome);std::cout<<genome<<std::endl;return0;}

Hints

Remember that the last element in genome is stored at position length(genome)-1.

Solution

Click more‚Ä¶ to see the solution.

#include<seqan/sequence.h>#include<seqan/basic.h>#include<seqan/stream.h>#include<seqan/file.h>#include<seqan/modifier.h>usingnamespaceseqan;DnagetRevCompl(Dnaconst&nucleotide){if(nucleotide=='A')return'T';if(nucleotide=='T')return'A';if(nucleotide=='C')return'G';return'C';}intmain(){DnaStringgenome="TATATACGCGCGAGTCGT";DnaStringrevComplGenome;//1.resize(revComplGenome,length(genome));//2.for(unsignedi=0;i<length(genome);++i)revComplGenome[length(genome)-1-i]=getRevCompl(genome[i]);//3.std::cout<<genome<<std::endl;std::cout<<revComplGenome<<std::endl;// And to check if your output is correct,// use the given SeqAn function reverseComplement(),// which modifies the sequence in-place:reverseComplement(genome);std::cout<<genome<<std::endl;return0;}

In this assignment, you will do some simple string building tasks, and write a simple alignment of the given reads and chromosomes.
Use the given code template to solve these subtasks:

Assume we have mapped the reads to the positions 7, 100, 172, and 272 in ‚Äėchr1‚Äô.
Store these positions in another string ‚ÄėalignPosList‚Äô.

Build another String bsChr1 as a copy of chr1, and exchange every ‚ÄėC‚Äô with a ‚ÄėT‚Äô, as in a bisulfite treated genome.

Print alignments of the reads and chr1 (or bschr1) using the function printAlign and the string alignPosList.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Function to print simple alignment between two sequences with the same lengthtemplate<typenameTText1,typenameTText2>voidprintAlign(TText1const&genomeFragment,TText2const&read){std::cout<<"Alignment "<<std::endl;std::cout<<" genome : "<<genomeFragment<<std::endl;std::cout<<" read : "<<read<<std::endl;}intmain(){// Build reads and genomesDnaStringchr1="TATAATATTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCTAGATGTGCAGCTCGATCGAATGCACGTGTGTGCGATCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATTTAG";// Build List containing all readstypedefString<DnaString>TDnaList;TDnaListreadList;resize(readList,4);readList[0]="TTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCT";readList[1]="TCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATT";readList[2]="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACA";readList[3]="CGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACA";// Append a second chromosome sequence fragment to chr1DnaStringchr2="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACACGTCTCTGTGTTCCGACGTGTGTCACGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACACATGCTGCTG";append(chr1,chr2);// Print readliststd::cout<<" \n Read list: "<<std::endl;for(unsignedi=0;i<length(readList);++i)std::cout<<readList[i]<<std::endl;// 1. Assume we have mapped the 4 reads to chr1 (and chr2) and now have the mapping start positions (no gaps).// Store the start position in a String alignPosList: 7, 100, 172, 272

// Your code snippet here for 1.+2.

// 3. Print alignments of the reads with chr1 (or bsChr1) sequence using the function printAlign// and the positions in alignPosList.// To do that, you have to create a copy of the fragment in chr1 (bsChr1) that is aligned to the read.std::cout<<" \n Print alignment: "<<std::endl;for(unsignedi=0;i<length(readList);++i){// Begin position beginPosition of a given alignment between the read and the genome

// Your code snippet here for 3.

// Genome fragmentDnaStringgenomeFragment;

// Your code snippet here for 3.

// Call of our function to print the simple alignmentprintAlign(genomeFragment,readList[i]);}return0;}

Hints

You have to create a copy of the fragment in chr1 (bsChr1) that is aligned to the read.

Solution

Click more‚Ä¶ to see the solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Function to print simple alignment between two sequences with the same lengthtemplate<typenameTText1,typenameTText2>voidprintAlign(TText1const&genomeFragment,TText2const&read){std::cout<<"Alignment "<<std::endl;std::cout<<" genome : "<<genomeFragment<<std::endl;std::cout<<" read : "<<read<<std::endl;}intmain(){// Build reads and genomesDnaStringchr1="TATAATATTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCTAGATGTGCAGCTCGATCGAATGCACGTGTGTGCGATCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATTTAG";// Build List containing all readstypedefString<DnaString>TDnaList;TDnaListreadList;resize(readList,4);readList[0]="TTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCT";readList[1]="TCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATT";readList[2]="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACA";readList[3]="CGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACA";// Append a second chromosome sequence fragment to chr1DnaStringchr2="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACACGTCTCTGTGTTCCGACGTGTGTCACGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACACATGCTGCTG";append(chr1,chr2);// Print readliststd::cout<<" \n Read list: "<<std::endl;for(unsignedi=0;i<length(readList);++i)std::cout<<readList[i]<<std::endl;// 1. Assume we have mapped the 4 reads to chr1 (and chr2) and now have the mapping start positions (no gaps).// Store the start position in a String alignPosList: 7, 100, 172, 272String<unsigned>alignPosList;resize(alignPosList,4);alignPosList[0]=7;alignPosList[1]=100;alignPosList[2]=172;alignPosList[3]=272;// 2. Bisulfite conversion// Assume chr1 is beeing bisulfate treated: Copy chr1 to a new genome bsChr1 and exchange every 'C' with a 'T'DnaStringbsChr1;assign(bsChr1,chr1);for(unsignedi=0;i<length(bsChr1);++i)if(bsChr1[i]=='C')bsChr1[i]='T';// 3. Print alignments of the reads with chr1 (or bsChr1) sequence using the function printAlign// and the positions in alignPosList.// To do that, you have to create a copy of the fragment in chr1 (bsChr1) that is aligned to the read.std::cout<<" \n Print alignment: "<<std::endl;for(unsignedi=0;i<length(readList);++i){// Begin position beginPosition of a given alignment between the read and the genomeunsignedbeginPosition=alignPosList[i];// Genome fragmentDnaStringgenomeFragment;// We have to create a copy of the corresponding fragment of the genome, where the read aligns tofor(unsignedj=0;j<length(readList[i]);++j)appendValue(genomeFragment,chr1[beginPosition+j]);// Call of our function to print the simple alignmentprintAlign(genomeFragment,readList[i]);}return0;}

Each comparison involves a scan of the two sequences for searching the first mismatch between the strings.
This could be costly if the two sequences share a long common prefix.
Suppose we want to branch in a program depending on whether a<b, a==b, or a>b.

In this case, although only one scan would be enough to decide what case is to be applied, each operator > and < performs a new comparison.
SeqAn offers the class Lexical to avoid unnecessary sequence scans.
Lexicals can store the result of a comparison, for example:

Move conversion.
If the source sequence is not needed any more after the conversion, it is always advisable to use move instead of assign.
The function move does not make a copy but can reuse the source sequence storage.
In some cases, move can also perform an in-place conversion.

In this assignment you will sort nucleotides.
Copy the code below. Adjust the code such that all nucleotides, which are lexicographically smaller than a Dna5 'G' are stored in a list lesser, while all nucleotides which are greater, should be stored in a list greater.
Print out the final lists.

#include<seqan/stream.h>#include<seqan/sequence.h>#include<seqan/file.h>usingnamespaceseqan;intmain(){String<Dna5>nucleotides="AGTCGTGNNANCT";String<Dna5>selected;// Append all elements of nucleotides, apart of Gs,// to the list selected.for(unsignedi=0;i<length(nucleotides);++i){appendValue(selected,nucleotides[i]);}std::cout<<"Selected nucleotides: "<<selected<<std::endl;return0;}

In this task you will compare whole sequences.
Reuse the code from above. Instead of a String<Dna5> we will now deal with a String<Dna5String>.
Build a string which contains the Dna5Strings ‚ÄúATATANGCGT‚ÄĚ, ‚ÄúAAGCATGANT‚ÄĚ and ‚ÄúTGAAANTGAC‚ÄĚ.
Now check for all elements of the container, if they are lexicographically smaller or bigger than the given subject sequence ‚ÄúGATGCATGAT‚ÄĚ and append them to a appropriate list.
Print out the final lists.

Hints

Try to avoid unnecessary sequence scans.

Solution

Click more‚Ä¶ to see the solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;intmain(){String<Dna5String>nucleotidesList;Dna5Stringstr1="ATATANGCGT";Dna5Stringstr2="AAGCATGANT";Dna5Stringstr3="TGAAANTGAC";resize(nucleotidesList,3);nucleotidesList[0]=str1;nucleotidesList[1]=str2;nucleotidesList[2]=str3;String<Dna5String>lesser;String<Dna5String>greater;Dna5Stringref="GATGCATGAT";// For each Dna5String of the String:for(unsignedi=0;i<length(nucleotidesList);++i){// Compare the Dna5String with the given reference string// The result of the comparison is stored in compLexical<>comp(nucleotidesList[i],ref);// The function isLess checks only the stored result// without comparing the sequences againif(isLess(comp))appendValue(lesser,nucleotidesList[i]);elseif(isGreater(comp))appendValue(greater,nucleotidesList[i]);}// Print the resultsstd::cout<<"Lesser sequences: "<<std::endl;for(unsignedi=0;i<length(lesser);++i){std::cout<<lesser[i]<<", ";}std::cout<<std::endl;std::cout<<"Greater sequences: "<<std::endl;for(unsignedi=0;i<length(greater);++i){std::cout<<greater[i]<<", ";}}

Very often you will be required to iterate over your string to either retrieve what‚Äôs stored in the string or to write something at a specific position.
For this purpose SeqAn provides Iterators for all container types.
The metafunction Iterator can be used to determine the appropriate iterator type for a given a container.

An iterator always points to one value of the container.
The operator operator* can be used to access this value by reference.
Functions like operator++(prefix) or operator‚Äď(prefix) can be used to move the iterator to other values within the container.

The functions begin and end, applied to a container, return iterators to the begin and to the end of the container.
Note that similar to C++ standard library iterators, the iterator returned by end does not point to the last value of the container but to the position behind the last one.
If the container is empty then end()==begin().

The following code prints out a sequence and demonstrates how to iterate over a string.

Some containers offer several kinds of iterators, which can be selected by an optional template parameter of the Iterator class.
For example, the tag Standard can be used to get an iterator type that resembles the C++ standard random access iterator.
For containers there is also a second variant available, the so called Rooted iterator.
The rooted iterator knows its container by pointing back to it.
This gives us a nice interface to access member functions of the underlying container while operating on a rooted iterator.
The construction of an iterator in SeqAn, e.g. for a Dna String, could look like the following:

Iterator<DnaString>::Typeit1;// A standard iteratorIterator<DnaString,Standard>::Typeit2;// Same as aboveIterator<DnaString,Rooted>::Typeit3;// A rooted iterator

Tip

The default iterator implementation is Standard.
Rooted iterators offer some convenience interfaces for the user.
They offer additional functions like container for determining the container on which the iterator works, and they simplify the interface for other functions like atEnd.
Moreover, rooted iterators may change the container‚Äôs length or capacity, which makes it possible to implement a more intuitive variant of a remove algorithm.

While rooted iterators can usually be converted into standard iterators, it is not always possible to convert standard iterators back into rooted iterators, since standard iterators may lack the information about the container they work on.
Therefore, many functions that return iterators like begin or end return rooted iterators instead of standard iterators; this way, they can be used to set both rooted and standard iterator variables.
Alternatively it is possible to specify the returned iterator type explicitly by passing the iterator kind as a tag argument, e.g. begin(str,Standard()).

Each sequence object has a capacity, i.e. the reserved space for this object.
The capacity can be set explicitly by functions such as reserve or resize.
It can also bet set implicitly by functions like append, assign, insert or replace, if the operation‚Äôs result exceeds the length of the target sequence.

If the current capacity of a sequence is exceeded by chaining the length, we say that the sequence overflows.
There are several overflow strategies that determine what actually happens when a string should be expanded beyond its capacity.
The user can specify this for a function call by additionally handing over a tag.
If no overflow strategy is specified, a default overflow strategy is selected depending on the type of the sequence.

Whenever the capacity is exceeded, the new capacity is chosen somewhat larger than currently needed.
This way, the number of capacity changes is limited in a way that resizing the sequence only takes amortized constant time.

No capacity check is performed, so the user has to ensure that the container‚Äôs capacity is large enough.

The next example illustrates how the different strategies could be used:

String<Dna>dnaSeq;// Sets the capacity of dnaSeq to 5.resize(dnaSeq,4,Exact());// Only "TATA" is assigned to dnaSeq, since dnaSeq is limited to 4.assign(dnaSeq,"TATAGGGG",Limit());std::cout<<dnaSeq<<std::endl;// Use the default expansion strategy.append(dnaSeq,"GCGCGC");std::cout<<dnaSeq<<std::endl;

Build a string of Dna (default specialization) and use the function appendValue to append a million times the nucleotide ‚ÄėA‚Äô.
Do it both using the overflow strategy Exact and Generous.
Measure the time for the two different strategies.

The user can specify the kind of string that should be used in an optional second template argument of String.
The default string implementation is Alloc String.

In most cases, the implementation Alloc String (the default when using a String<T>) is the best choice.
Exceptions are when you want to process extremely large strings that are a bit larger than the available memory (consider Alloc String) or much larger so most of them are stored on the hard disk and only parts of them are loaded in main memory (consider External String).
The following list describes in detail the different specializations:

// Most of the string is stored on the disk.String<Dna,External<>>myLargeGenome;

Tip

String Simplify Memory Management

One advantage of using Strings is that the user does not need to reserve memory manually with new and does not need delete to free memory.
Instead, those operations are automatically handled by the String class.

The following section will introduce you into the Segment class of SeqAn.

Segments are contiguous subsequences that represent parts of other sequences.
Therefore, their functionality is similar to the String functionality.
In SeqAn, there are three kinds of segments: InfixSegment, PrefixSegment, and SuffixSegment.
The metafunctions Infix, Prefix, and Suffix, respectively, return the appropriate segment data type for a given sequence type.

For prefixes, we use the function prefix to build the prefix.
The first parameter is the sequence we build the prefix from, the second the excluding end position.
For infixes, we have to provide both the including start and the excluding end position.
For suffixes, the second parameter of the function denotes the including starting position of the suffix:

Segments store a pointer on the underlying sequence object, the host, and an start and/or end position, depending on the type of segment.
The segment is not a copy of the sequence segment.

Warning

Please note that it is not possible anymore to change the underlying sequence by changing the segment.
If you want to change the host sequence, you have to explicitly modify this.
If you want to modify only the segment, you have to explicitly make a copy of the string.

In this task you will use a segment to pass over an infix of a given sequence to a function without copying the corresponding fragment.
Use the code given below.
Lets assume that we have given a genome and a read sequence as well as the begin position of a given alignment.
In the main function a fragment of the Dna5String genome is copied and passed together with the Dna5String read to a print function.
Adjust the code to use an infix of the genome, instead of copying the corresponding fragment.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Function to print simple alignment between two sequences with the same length// .. for two sequences of different typestemplate<typenameTText1,typenameTText2>voidprintAlign(TText1const&genomeFragment,TText2const&read){std::cout<<"Alignment "<<std::endl;std::cout<<" genome : ";std::cout<<genomeFragment<<std::endl;std::cout<<" read : ";std::cout<<read<<std::endl;}intmain(){// We have given a genome sequenceDna5Stringgenome="ATGGTTTCAACGTAATGCTGAACATGTCGCGT";// A read sequenceDna5Stringread="TGGTNTCA";// And the begin position of a given alignment between the read and the genomeunsignedbeginPosition=1;

// We have to create a copy of the corresponding fragment of the genome, where the read aligns to// Change this piece of code using an infix of the genomeDna5StringgenomeFragment;for(unsignedi=0;i<length(read);++i){appendValue(genomeFragment,genome[beginPosition+i]);}

// Call of our function to print the simple alignmentprintAlign(genomeFragment,read);return0;}

Solution

Click more‚Ä¶ to see the solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Function to print simple alignment between two sequences with the same length// .. for two sequences of different typestemplate<typenameTText1,typenameTText2>voidprintAlign(TText1const&genomeFragment,TText2const&read){std::cout<<"Alignment "<<std::endl;std::cout<<" genome : ";std::cout<<genomeFragment<<std::endl;std::cout<<" read : ";std::cout<<read<<std::endl;}intmain(){// We have given a genome sequenceDna5Stringgenome="ATGGTTTCAACGTAATGCTGAACATGTCGCGT";// A read sequenceDna5Stringread="TGGTNTCA";// And the begin position of a given alignment between the read and the genomeunsignedbeginPosition=1;// Create Infix of type Dna5String and get the corresponding infix sequence of genomeInfix<Dna5String>::TypegenomeFragment=infix(genome,beginPosition,beginPosition+length(read));// Call of our function to print the simple alignmentprintAlign(genomeFragment,read);return0;}

Take the solution from the workshop assignment above and change it to use Segments for building the genome fragment.

Hints

Note that because printAlign uses templates, you don‚Äôt have to change the function even though the type of genomeFragment is different.

Solution

Click more‚Ä¶ to see the solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Function to print simple alignment between two sequences with the same lengthtemplate<typenameTText1,typenameTText2>voidprintAlign(TText1const&genomeFragment,TText2const&read){std::cout<<"Alignment "<<std::endl;std::cout<<" genome : "<<genomeFragment<<std::endl;std::cout<<" read : "<<read<<std::endl;}intmain(int,charconst**){// Build reads and genomesDnaStringchr1="TATAATATTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCTAGATGTGCAGCTCGATCGAATGCACGTGTGTGCGATCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATTTAG";// Build List containing all readstypedefString<DnaString>TDnaList;TDnaListreadList;resize(readList,4);readList[0]="TTGCTATCGCGATATCGCTAGCTAGCTACGGATTATGCGCTCTGCGATATATCGCGCT";readList[1]="TCGATTAGCGTCGATCATCGATCTATATTAGCGCGCGGTATCGGACGATCATATTAGCGGTCTAGCATT";readList[2]="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACA";readList[3]="CGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACA";// Append a second chromosome sequence fragment to chr1DnaStringchr2="AGCCTGCGTACGTTGCAGTGCGTGCGTAGACTGTTGCAAGCCGGGGGTTCATGTGCGCTGAAGCACACATGCACACGTCTCTGTGTTCCGACGTGTGTCACGTGCACTGCTGACGTCGTGGTTGTCACATCGTCGTGCGTGCGTACTGCTGCTGACACATGCTGCTG";append(chr1,chr2);// Print readliststd::cout<<" \n Read list: "<<std::endl;for(unsignedi=0;i<length(readList);++i)std::cout<<readList[i]<<std::endl;// Assume we have mapped the 4 reads to chr1 (and chr2) and now have the mapping start positions (no gaps).// Store the start position in a String alignPosList: 7, 100, 172, 272String<unsigned>alignPosList;resize(alignPosList,4);alignPosList[0]=7;alignPosList[1]=100;alignPosList[2]=172;alignPosList[3]=272;// Optional// Bisulfite conversion// Assume chr1 is beeing bisulfate treated: Copy chr1 to a new genome bsChr1 and exchange every 'C' with a 'T'DnaStringbsChr1;assign(bsChr1,chr1);for(unsignedi=0;i<length(bsChr1);++i)if(bsChr1[i]=='C')bsChr1[i]='T';// Print alignments using Segment: Do the same as above, but instead of using a for loop to build the fragment,// use the Segment class to build an infix of bsChr1.// Note: Because printAlign uses templates, we don't have to change the function even though the type of// genomeFragment is different.std::cout<<" \n Print alignment using Segment: "<<std::endl;for(unsignedi=0;i<length(readList);++i){// Begin and end position of a given alignment between the read and the genomeunsignedbeginPosition=alignPosList[i];unsignedendPosition=beginPosition+length(readList[i]);// Build infixInfix<DnaString>::TypegenomeFragment=infix(chr1,beginPosition,endPosition);// Call of our function to print the simple alignmentprintAlign(genomeFragment,readList[i]);}return0;}

A set of sequences can either be stored in a sequence of sequences, for example in a String<String<char>>, or in a StringSet.
This tutorial will introduce you to the SeqAn class StringSet, its background and how to use it.

One advantage of using StringSet is that it supports the function concat that returns a concatenator of all sequences in the string set.
A concatenator is an object that represents the concatenation of a set of strings.
This way, it is possible to build up index data structures for multiple sequences by using the same construction methods as for single sequences.

StringSet<DnaString>ownerSet;StringSet<DnaString,Owner<>>ownerSet2;// same as aboveStringSet<DnaString,Dependent<>>dependentSet;

The specialization ConcatDirect StringSet already stores the sequences in a concatenation.
The concatenators for all other specializations of StringSet are virtual sequences, that means their interface simulates a concatenation of the sequences, but they do not literally concatenate the sequences into a single sequence.
Hence, the sequences do not need to be copied when a concatenator is created.

One string can be an element of several Dependent StringSets.
Typical tasks are, e.g., to find a specific string in a string set, or to test whether the strings in two string sets are the same.
Therefore a mechanism to identify the strings in the string set is needed, and, for performance reasons, this identification should not involve string comparisons.
SeqAn solves this problem by introducing ids, which are by default unsignedint values.

The sequences are stored as parts of a long string.
Since the sequences are already concatenated, concat just needs to return this string.
The string set also stores lengths and starting positions of the strings.
Inserting new strings into the set or removing strings from the set is more expensive than for the default OwnerStringSet specialization, since this involves moving all subsequent sequences in memory.

Specialization Owner<JournaledSet>

The sequences are stored as Journaled Strings to a common reference sequence,
that is also stored within the container.
When adding a new String to the set, it needs to be joined to this set of sequences which are all based on the
common reference sequence.
This way one can hold a large collection of similar sequences efficiently in memory.

Specialization Dependent<Tight>

This specialization stores sequence pointers consecutively in an array.
Another array stores an id value for each sequence.
That means that accessing given an id needs a search through the id array.

Warning

The Dependent-Tight StringSet is deprecated and will likely be removed
within the SeqAn-2.x lifecycle.

Specialization Dependent<Generous>

The sequence pointers are stored in an array at the position of their ids.
If a specific id is not present, the array stores a zero at this position.
The advantage of this specialization is that accessing the sequence given its id is very fast.
On the other hand, accessing a sequence given its position i can be expensive, since this means we have to find the i-th non-zero value in the array of sequence pointers.
The space requirements of a string set object depends on the largest id rather than the number of sequences stored in the set.
This could be inefficient for string sets that store a small subset out of a large number of sequences.

This section will give you a short overview of the functionality of the class StringSet.

There are two ways for accessing the sequences in a string set: (1) the function operator[] returns a reference to the sequence at a specific position within the sequence of sequences, and (2) valueById accesses a sequence given its id.
We can retrieve the id of a sequence in a StringSet with the function positionToId.

// Let's create a string set of type dependent to represent strings,// which are stored in the StringSet of type OwnerStringSet<DnaString,Dependent<Tight>>depSet;// We assign the first two strings of the owner string set to the dependent StringSet,// but in a reverse orderassignValueById(depSet,stringSet,id1);assignValueById(depSet,stringSet,id0);std::cout<<"Dependent: "<<'\n';// (1) Access by positionstd::cout<<"Pos 0: "<<value(depSet,0)<<'\n';// (2) Access by idstd::cout<<"Id 0: "<<valueById(depSet,id0)<<'\n';

Dependent:Pos 0: CGCGId 0: TATA

With the function positionToId we can show that, in this case, the position and the id of a string are different.

If we want to iterate over the contained Strings as well, as if the StringSet would be one sequence, we can use the function concat to get the concatenation of all sequences.
Therefore we first use the metafunction Concatenator to receive the type of the concatenation.
Then, we can simply build an iterator for this type and iterate over the concatenation of all strings.

Build a string set with default specialization and which contains the strings "AAA", "CCC", "GGG" and "TTT".
After that print the length of the string set and use a simple for-loop to print all elements of the strings set.

Solution

Click more‚Ä¶ to see the solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;intmain(){// Build stringsDnaStringstr0="AAA";DnaStringstr1="CCC";DnaStringstr2="GGG";DnaStringstr3="TTT";// Build string set and append stringsStringSet<DnaString>stringSet;appendValue(stringSet,str0);appendValue(stringSet,str1);appendValue(stringSet,str2);appendValue(stringSet,str3);// Print the length of the string setstd::cout<<length(stringSet)<<std::endl;// Print all elementsfor(unsignedi=0;i<length(stringSet);++i){std::cout<<stringSet[i]<<std::endl;}return0;}

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>usingnamespaceseqan;// Check whether the string set contains the string with the given id,// without comparing the actual sequencestemplate<typenameTStringSet,typenameTId>boolisElement(TStringSet&stringSet1,TId&id){for(unsignedi=0;i<length(stringSet1);++i){// Get the id of the element at position iif(positionToId(stringSet1,i)==id)returntrue;}returnfalse;}intmain(){// Build stringsDnaStringstr0="TATA";DnaStringstr1="CGCG";DnaStringstr2="TTAAGGCC";DnaStringstr3="ATGC";DnaStringstr4="AGTGTCA";// Build owner string set and append stringsStringSet<DnaString>stringSetOw;appendValue(stringSetOw,str0);appendValue(stringSetOw,str1);appendValue(stringSetOw,str2);appendValue(stringSetOw,str3);appendValue(stringSetOw,str4);// Get corresponding ids for positionsunsignedid0=positionToId(stringSetOw,0);unsignedid1=positionToId(stringSetOw,1);unsignedid2=positionToId(stringSetOw,2);unsignedid3=positionToId(stringSetOw,3);// Build dependent string set and assigns strings by idStringSet<DnaString,Dependent<Generous>>stringSetDep;assignValueById(stringSetDep,stringSetOw,id0);assignValueById(stringSetDep,stringSetOw,id1);assignValueById(stringSetDep,stringSetOw,id3);// Call function to check if a string is contained and output resultstd::cout<<"Does the string set contain the string with the id 'id3'? "<<isElement(stringSetDep,id3)<<std::endl;std::cout<<"Does the string set contain the string with the id 'id2'? "<<isElement(stringSetDep,id2)<<std::endl;return0;}

This tutorial will describe the different alphabets used in SeqAn, or in other words, you will learn about the contained types of a SeqAn String.
To continue with the other tutorials, it would be enough to know, that in SeqAn several standard alphabets are already predefined, e.g. Dna, Dna5, Rna, Rna5, Iupac, AminoAcid.

Any type that provides a default constructor, a copy constructor and an assignment operator can be used as the alphabet / contained type of a String (see also the tutorial Sequences).
This includes the C++ POD types, e.g. char, int, double etc.
In addition you can use more complex types like String as the contained type of strings, e.g. String<String<char>>.

SeqAn also provides the following types that are useful in bioinformatics.
Each of them is a specialization of the class SimpleType.

In SeqAn, alphabets are value types that can take a limited number of values and which hence can be mapped to a range of natural numbers.
We can retrieve the number of different values of an alphabet, the alphabet size, by the metafunction ValueSize.

Another useful metafunction called BitsPerValue can be used to determine the number of bits needed to store a value of a given alphabet.

unsignedbits=BitsPerValue<TAlphabet>::VALUE;std::cout<<"Number of bits needed to store a value of type Dna: "<<bits<<'\n';

Number of bits needed to store a value of type Dna: 2

The order of a character in the alphabet (i.e. its corresponding natural number) can be retrieved by calling the function ordValue.
See each specialization‚Äôs documentation for the ordering of the alphabet‚Äôs values.

The return value of the ordValue function is determined by the metafunction ValueSize.
ValueSize returns the type which uses the least amount of memory while being able to represent all possible values.
E.g. ValueSize of Dna returns an _uint8 which is able to represent 256 different characters.
However, note that std::cout has no visible symbol for printing all values on the screen, hence a cast to unsigned might be necessary.

In this task you will learn how to access all the letters of an alphabet.
Use the piece of code from below and adjust the function showAllLettersOfMyAlphabet() to go through all the characters of the current alphabet and print them.

#include<seqan/sequence.h>#include<seqan/basic.h>#include<iostream>usingnamespaceseqan;// We want to define a function, which takes// the alphabet type as an argumenttemplate<typenameTAlphabet>voidshowAllLettersOfMyAlphabet(TAlphabetconst&){// ...}intmain(){showAllLettersOfMyAlphabet(AminoAcid());showAllLettersOfMyAlphabet(Dna());showAllLettersOfMyAlphabet(Dna5());return0;}

#include<seqan/sequence.h>#include<seqan/basic.h>#include<iostream>usingnamespaceseqan;// We define a function which takes// the alphabet type as an argumenttemplate<typenameTAlphabet>voidshowAllLettersOfMyAlphabet(TAlphabetconst&){typedeftypenameValueSize<TAlphabet>::TypeTSize;// We need to determine the alphabet size// using the metafunction ValueSizeTSizealphSize=ValueSize<TAlphabet>::VALUE;// We iterate over all characters of the alphabet// and output themfor(TSizei=0;i<alphSize;++i)std::cout<<static_cast<unsigned>(i)<<','<<TAlphabet(i)<<" ";std::cout<<std::endl;}intmain(){showAllLettersOfMyAlphabet(AminoAcid());showAllLettersOfMyAlphabet(Dna());showAllLettersOfMyAlphabet(Dna5());return0;}

Tasks, such as computing an alignment, searching for patterns online or indexing a genome or protein database are required in many bioinformatics applications.
For these tasks, SeqAn provides fundamental data structures to work efficiently with biological sequences.
SeqAn implements special alphabets, such as DNA, RNA, AminoAcid, Iupac and more.
The alphabets available in SeqAn can be reviewed here.
You will also find some more information about using the alphabet types.

Besides the alphabets SeqAn also implements a string class.
SeqAn strings are generic containers in which characters of any alphabet are stored continuously in memory.
The default string class implementation is equivalent to the STL vector class.
However, the memory mangement of the SeqAn string class is optimized for working with SeqAn‚Äôs alphabets.
Apart of the default string class implementation SeqAn provides many useful specializations of this class, which are very useful in the bioinformatics context.
The tutorial about Strings and Segments gives you a more detailed overview over the string class and it‚Äôs functionality.

Another generic container data structure used very often is the string set.
The string set is a special container to represent a collection of strings, which in addition provides many helpful functions to make the work even easier and more efficient.
The StringSet tutorial introduces you to this class and how to work with it.

An FMIndex is a space-efficient and fast index structure that can be used to search arbitrary patterns.
In contrast to other indices like the IndexSa, the FMIndex can only be traversed as a prefix trie.
Patterns cannot be searched from left to right but from right to left.
This means that you either have to reverse the text before constructing the index or reverse the patterns before searching.
Also hits of a pattern do not refer to the starting position of the pattern but to the ending position.
This has to be taken into account when retrieving hits.

The FMIndex is based on the Burrow-Wheeler-Transform (BWT) which uses a rank dictionary.
We currently have three different implementations to choose from with different running time and space trade-offs (the space consumption only refers to the rank dictionary, not the entire index).
\(n\) is the length of the text, \(\sigma\) the alphabet size.
The running time refers to a single character search.

Generally speaking, WaveletTrees are recommended for larger alphabets, while Levels offer a faster runtime at the expense of a higher space consumption.
LevelsRDConfigs are in practice noticeably faster than WaveletTrees since they only perform \(\mathcal{O}(\log \sigma)\) logical operations instead of going down a tree of height \(\mathcal{O}(\log \sigma)\) resulting in more cache misses.
LevelsPrefixRDConfig are recommended for smaller alphabets and bidirectional support on FM indices.
The running time is constant for both backward and forward searches.

All three rank dictionaries are based on bit vectors with constant-time rank support.
We skip at this point details on how rank queries work and focus on the parameters and their effects.
To reduce the space consumption there are three parameters that you can adjust:

First of all TSize should be the smallest data type that can store the length of the indexed text (e.g. uint8_t, uint16_t, uint32_t, uint64_t).
For a StringSet the length is defined as the the sum of lengths of strings plus the number of strings.

The bit vector is clustered into blocks that store the precomputed rank up to every block.
To reduce the space of these blocks, one can add an additional level of blocks on top which group blocks together to a superblock.
We support up to 3 levels.
The effect on the running time is minimal since for every additional level only one more array lookup is conducted.
In practice two levels are in most cases faster than one level due to smaller tables and thus less cache misses.

To reduce the space consumption even further or to improve the running time, one can change the size of a block.
Each block contains \(64 \cdot WPB\) (WORDS_PER_BLOCK) bits and its rank can thus be computed on a 64 bit machine by \(WPB\) popcount operations.
By reducing (increasing) the block size, the running time (space consumption) can be improved noticeably.
By default, WPB is set to 1.

For texts with less than \(2^{32}\) characters, two levels and a fast query time (only one popcount operation per rank query), the FMIndex can be configured using

The bidirectional FMIndex can be used to search a pattern into both directions, i.e. a string can be extend by a character to the left or right in an arbitrary manner.
For information on how to search in a bidirectional FMIndex, please check the Tutorial on Index Iterators.

A q-gram index can be used to efficiently retrieve all occurrences of a certain q-gram in the text.
It consists of various tables, called fibres (see Accessing Index Fibres), to retrieve q-gram positions, q-gram counts, etc.
However, it has no support for suffix tree iterators.
A q-gram index must be specialized with a Shape type.
A Shape defines q, the number of characters in a q-gram and possibly gaps between these characters.
There are different specializations of Shape available:

Each shape evaluates a gapped or ungapped sequence of q characters to a hash value by the Functions hash, hashNext, etc.
For example, the shape 1101 represents a 3-gram with one gap of length 1.
This shape overlayed with the Dna text "GATTACA" at the third position corresponds to "TT-C".
The function hash converts this 3-gram into \(61 = (\mathbf{3} \cdot 4 + \mathbf{3}) \cdot 4 + 1\).
4 is the alphabet size in this example (see ValueSize).

With hash and hashNext, we can compute the hash values of arbitrary / adjacent q-grams and a loop that outputs the hash values of all overlapping ungapped 3-grams could look as follows:

Note that the shape not only stores the length and gaps of a q-gram shape but also stores the hash value returned by the last hash/hashNext call.
This hash value can be retrieved by calling value on the shape.
However, one drawback of the example loop above is that the first hash value must be computed with hash while the hash values of the following overlapping q-grams can more efficiently be computed by hashNext.
This complicates the structure of algorithms that need to iterate all hash values, as they have to handle this first hash differently.
As a remedy, the hashInit function can be used first and then hashNext on the first and all following text positions in the same way:

The q-gram index offers different functions to search or count occurrences of q-grams in an indexed text, see getOccurrences, countOccurrences.
A q-gram index over a StringSet stores occurrence positions in the same way and in the same fibre (FibreSA) as the ESA index.
If only the number of q-grams per sequence are needed the QGramCounts and QGramCountsDir fibres can be used.
They store pairs (seqNo,count), count>0, for each q-gram that occurs counts times in sequence number seqNo.

To efficiently retrieve all occurrence positions or all pairs (seqNo,count) for a given q-gram, these positions or pairs are stored in contiguous blocks (in QGramSA, QGramCounts fibres), called buckets.
The begin position of bucket i is stored in directory fibres (QGramDir, QGramCountsDir) at position i, the end position is the begin positions of the bucket i+1.
The default implementation of the IndexQGram index maps q-gram hash values 1-to-1 to bucket numbers.
For large q or large alphabets the Open Addressing QGram Index can be more appropriate as its directories are additionally bound by the text length.
This is realized by a non-trivial mapping from q-gram hashes to bucket numbers that requires an additional fibre (QGramBucketMap).

We want to construct the q-gram index of the string "CATGATTACATA" and output the occurrences of the ungapped 3-gram "CAT".
As 3 is fixed at compile-time and the shape has no gaps we can use an UngappedShape which is the first template argument of IndexQGram, the second template argument of Index.
Next we create the string "CATGATTACATA" and specialize the first index template argument with the type of this string.
The string can be given to the index constructor.

To get all occurrences of a q-gram, we first have to hash it with a shape of the same type as the index shape (we can even use the index shape returned by indexShape).
The hash value returned by hash or hashNext is also stored in the shape and is used by the function getOccurrences to retrieve all occurrences of our 3-gram.

Write a program that outputs all occurrences of the gapped q-gram ‚ÄúAT-A‚ÄĚ in ‚ÄúCATGATTACATA‚ÄĚ.

Solution

Before we can create a DnaString index of ‚ÄúCATGATTACATA‚ÄĚ, we have to choose an appropriate Shape.
Because our shape 1101 is known at compile-time and contains only one gap we could choose OneGappedShape, GappedShape, or GenericShape (see the commented-out code).
Although the GenericShape could be used for every possible shape, it is a good idea to choose a Shape with restrictions as its hash functions are more efficient in general.

Please note that the Shape object that corresponds to the IndexQGram index is empty initially and has to be set by stringToShape or resize.
This initialization is not necessary for Shape that are defined at compile-time, i.e. UngappedShape and GappedShape.
To search for ‚ÄúAT-A‚ÄĚ we first have to hash it with the index shape or any other Shape with the same bitmap.
The we can use getOccurrences to output all matches.

Instead of length(getOccurrences(...)) we could have used countOccurrences.
But beware that countOccurrences requires only the QGram_Dir fibre, whereas getOccurrences requires both QGram_Dir and QGram_SA, see Accessing Index Fibres.
Because QGram_SA can be much more efficiently constructed during the construction of QGram_Dir, QGram_Dir would be constructed twice.

Create and output a matrix M where M(i,j) is the number of common ungapped 5-grams between sequence i and sequence j for 3 random Dna sequences, each not longer than 200 characters.
Optional: Run the matrix calculation twice, once for an IndexQGram and once for an Open Addressing QGram Index and output the directory sizes (QGram_Dir, QGram_CountsDir fibre).

Hint

A common q-gram that occurs \(a\) times in one and \(b\) times in the other sequence counts for \(\min(a,b)\).

Solution

For generating random numbers we use the std::mt19937.
The random numbers returned by the random number engine are arbitrary unsignedint values which we downscale to values between 0 and 3 and convert into Dna characters.
The 3 generated strings are of random length and appended to a StringSet.
The main algorithm is encapsulated in a template function qgramCounting to easily switch between the two IndexQGram specializations.

intmain(){// for the sake of reproducibilitystd::mt19937rng;// create StringSet of 3 random sequencesStringSet<DnaString>stringSet;reserve(stringSet,3);for(intseqNo=0;seqNo<3;++seqNo){DnaStringtmp;intlen=rng()%100+10;for(inti=0;i<len;++i)appendValue(tmp,Dna(rng()%4));appendValue(stringSet,tmp);std::cout<<">Seq"<<seqNo<<std::endl<<tmp<<std::endl;}qgramCounting(stringSet,IndexQGram<UngappedShape<5>>());qgramCounting(stringSet,IndexQGram<UngappedShape<5>,OpenAddressing>());return0;}

The main function expects the StringSet and the Index specialization as a tag.
First, we define lots of types we need to iterate and access the fibres directly.
We then notify the index about the fibres we require.
For storing the common q-grams we use a 2-dimensional Matrix object whose lengths have to be set with setLength for each dimension.
The matrix is initialized with zeros by resize.

The main part of the function iterates over the CountsDir fibre.
Each entry in this directory represents a q-gram bucket, a contiguous interval in the Counts fibre storing for every sequence the q-gram occurs in the number of occurrences in pairs (seqNo,count).
The interval begin of each bucket is stored in the directory and the interval end is the begin of the next bucket.
So the inner loops iterate over all non-empty buckets and two pairs (seqNo1,count1) and (seqNo2,count2) indicate that seqNo1 and seqNo2 have a common q-gram.
At the end the matrix can simply be output by shifting it to the cout stream.

// for each bucket count common q-grams for each sequence pairTDirValuebucketBegin=*itCountsDir;for(++itCountsDir;itCountsDir!=itCountsDirEnd;++itCountsDir){TDirValuebucketEnd=*itCountsDir;// q-gram must occur in at least 2 different sequencesif(bucketBegin!=bucketEnd){TIterCountsitA=itCountsBegin+bucketBegin;TIterCountsitEnd=itCountsBegin+bucketEnd;for(;itA!=itEnd;++itA)for(TIterCountsitB=itA;itB!=itEnd;++itB)distMat((*itA).i1,(*itB).i1)+=_min((*itA).i2,(*itB).i2);}bucketBegin=bucketEnd;}std::cout<<std::endl<<"Common 5-mers for Seq_i, Seq_j"<<std::endl;std::cout<<distMat;}

Indices in SeqAn are substring indices, meaning that they allow efficient pattern queries in strings or sets of strings.
In contrast to, e.g., online-search algorithms that search through the text in \(\mathcal{O}(n)\), substring indices find a pattern in sublinear time \(o(n)\).

We will now show how we can create the different indices in SeqAn before we show how they are used for pattern search.

All the mentioned indices belong to the generic Index class.
A SeqAn index needs two pieces of information: the type of the String or StringSet to be indexed and the index specialization, such as IndexEsa or FMIndex.

Important

Indices based on suffix arrays (also including the FM index) are built using secondary memory.
When building large indices, it is therefore possible to run out of disk space (in which case an exception will be
thrown).
To circumvent this, the directory used for temporary storage can be changed by specifying the TMPDIR environment variable (on UNIX)
respectively TEMP environment variable (on Windows):

#exportTMPDIR=/somewhere/else/with/more/space

# SET TEMP=C:\somewhere\else\with\more\space

The following code snippet creates an enhanced suffix array index of a string of type Dna5.

#include<seqan/sequence.h>#include<seqan/index.h>usingnamespaceseqan;intmain(){String<char>text="This is the first example";Index<String<char>,FMIndex<>>index(text);return0;}

Solution

#include<seqan/sequence.h>#include<seqan/index.h>usingnamespaceseqan;intmain(){// One possible solution to the first sub assignmentString<Dna>text="ACGTTTGACAGCT";Index<String<Dna>,IndexEsa<>>index(text);// One possible solution to the second sub assignmentStringSet<String<Dna>>stringSet;appendValue(stringSet,"ACGTCATCAT");appendValue(stringSet,"ACTTTG");appendValue(stringSet,"CACCCCCCTATTT");Index<StringSet<String<Dna>>,IndexEsa<>>indexSet(stringSet);return0;}

You might have noticed that we only applied the FMIndex and IndexEsa in the examples.
The reason for this is that even though everything stated so far is true for the other indices as well, IndexWotd and IndexDfi are more useful when used with iterators as explained in the tutorial Index Iterators and the IndexQGram uses Shapes which is also explained in another tutorial.

One last remark is necessary.

Important

If you search for two different patterns with the same Finder object, you have to call the clear function of the finder between the search for the two patterns.
Otherwise the behavior is undefined.

The previous sections already described how an index of a set of strings can be instantiated.
A character position of a StringSet can be one of the following:

A local position (default), i.e. a Pair (seqNo, seqOfs) where seqNo identifies the string within the StringSet and the seqOfs identifies the position within this string.

A global position, i.e. a single integer value between 0 and the sum of string lengths minus 1.
This integer is the position in the gapless concatenation of all strings in the StringSet to a single string.

For indices, the meta-function SAValue determines, which position type (local or global) will be used for internal index tables (suffix array, q-gram array) and what type of position is returned by functions like position of a Finder.
SAValue returns a Pair (local position) by default, but could be specialized to return an integer type (global position) for some applications.
If you want to write algorithms for both variants you should use the functions posLocalize, posGlobalize, getSeqNo, and getSeqOffset.

If you have built your q-gram index with variable shapes (i.e. SimpleShapeGenericShape), you have to keep in mind that q or the shape is not stored or loaded.
This must be done manually and directly before or after loading with resize or stringToShape.

A newly instantiated index is initially empty.
If you assign a text to be indexed, solely the text fibre is set.
All other fibres are empty and created on demand.
Normally, a full created index should be saved to disk.
Therefore, you have to create the required fibres explicitly by hand.

Indexes based on external strings, e.g. Index<String<Dna,External<>>,IndexEsa<>> or Index<String<Dna,MMap<>>,IndexEsa<>> cannot be saved, as they are persistent implicitly.
The first thing after instantiating such an index should be associating it to a file with:

The file association implies that any change on the index, e.g. fibre construction, is synchronized to disk.
When instantiating and associating the index the next time, the index contains its previous state and all yet to be constructed fibres.

All Indices in SeqAn are capable of indexing Strings or StringSets of arbitrary sizes, i.e. up to 2^64 characters.
This always comes at a cost in terms of memory consumption, as any Index has to represent 64 bit positions in the underlying text.
However, in many practical instances, the text to be indexed is shorter, e.g. it does not exceed 4.29 billion (2^32) characters.
In this case, one can reduce the memory consumption of an Index by changing its internal data types, with no drawback concerning running time.

All Indices in SeqAn internally use the FibreSA, i.e. some sort of suffix array.
For Strings, each suffix array entry consumes 64 bit of memory per default, where 32 bit would be sufficient if the text size is appropriate.
In order to change the size type of the suffix array entry we simply have to overload the metafunction SAValue.

template<>structSAValue<String<Dna>>{typedefunsignedType;};

If your text is a StringSet, then SAValue will return a Pair that can be overloaded in the same way.

In the next example we want to use the TopDown Iterator to efficiently search a text for exact matches of a pattern.
We therefore want to use goDown which has an overload to go down an edge beginning with a specific character.

Important

The following examples show how to iterate IndexSa, IndexEsa, IndexWotd or IndexDfi, i.e. Index specializations representing suffix trees.
The result of the iteration will look different on Index specializations representing tries, e.g. FMIndex or IndexQGram.
Indeed, the topology of an Index changes depending on the chosen tree or trie specialization.
Note that any suffix tree edge can be labeled by more than one character, whereas any trie edge is always labeled by exactly one character.

First we create an index of the text "Howmuchwoodwouldawoodchuckchuck?"

intmain(){typedefIndex<CharString>TIndex;TIndexindex("How much wood would a woodchuck chuck?");Iterator<TIndex,TopDown<>>::Typeit(index);CharStringpattern="wood";while(repLength(it)<length(pattern)){// go down edge starting with the next pattern characterif(!goDown(it,pattern[repLength(it)]))return0;unsignedendPos=std::min((unsigned)repLength(it),(unsigned)length(pattern));// compare remaining edge characters with patternstd::cout<<representative(it)<<std::endl;if(infix(representative(it),parentRepLength(it)+1,endPos)!=infix(pattern,parentRepLength(it)+1,endPos))return0;}// if we get here the pattern was found// output match positionsfor(unsignedi=0;i<length(getOccurrences(it));++i)std::cout<<getOccurrences(it)[i]<<std::endl;return0;}

Afterwards we create the TopDown Iterator using the metafunction Iterator, which expects two arguments, the type of the container to be iterated and a specialization tag (see the VSTree Iterator hierarchy and the Iteration Tutorial for more details).

Iterator<TIndex,TopDown<>>::Typeit(index);

The main search can then be implemented using the functions repLength and representative.
Since goDown might cover more than one character (when traversing trees) it is necessary to compare parts of the pattern against the representative of the iterator.
The search can now be implemented as follows.
The algorithm descends the suffix tree along edges beginning with the corresponding pattern character.
In each step the unseen edge characters have to be verified.

CharStringpattern="wood";while(repLength(it)<length(pattern)){// go down edge starting with the next pattern characterif(!goDown(it,pattern[repLength(it)]))return0;unsignedendPos=std::min((unsigned)repLength(it),(unsigned)length(pattern));// compare remaining edge characters with patternstd::cout<<representative(it)<<std::endl;if(infix(representative(it),parentRepLength(it)+1,endPos)!=infix(pattern,parentRepLength(it)+1,endPos))return0;}

If all pattern characters could successfully be compared we end in the topmost node who‚Äôs leaves point to text positions starting with the pattern.
Thus, the suffixes represented by this node are the occurrences of our pattern and can be retrieved with getOccurrences.

// if we get here the pattern was found// output match positionsfor(unsignedi=0;i<length(getOccurrences(it));++i)std::cout<<getOccurrences(it)[i]<<std::endl;return0;}

Program output:

wwowood922

Alternatively, we could have used goDown to go down the path of the entire pattern instead of a single characters:

When implementing recursive algorithms such as an approximate search using backtracking, we recommend the use of the TopDownIterator without history.
By passing the iterator by value, the history is stored implicitly on the call stack.

Copy the code into a demo program and replace the text with a string set containing the strings "Howmuch",``‚ÄĚwood would‚ÄĚ`` and "awoodchuckchuck?".

Solution

#include<iostream>#include<seqan/index.h>usingnamespaceseqan;intmain(){StringSet<String<char>>text;appendValue(text,"How much");appendValue(text," wood would");appendValue(text," a woodchuck chuck?");typedefIndex<StringSet<String<char>>>TIndex;TIndexindex(text);Iterator<TIndex,TopDown<>>::Typeit(index);CharStringpattern="wood";while(repLength(it)<length(pattern)){// go down edge starting with the next pattern characterif(!goDown(it,pattern[repLength(it)]))return0;unsignedendPos=_min(repLength(it),length(pattern));// compare remaining edge characters with patternstd::cout<<representative(it)<<std::endl;if(infix(representative(it),parentRepLength(it)+1,endPos)!=infix(pattern,parentRepLength(it)+1,endPos))return0;}// if we get here the pattern was found// output match positionsfor(unsignedi=0;i<length(getOccurrences(it));++i)std::cout<<getOccurrences(it)[i]<<std::endl;return0;}

The difference is the format of the positions of the found occurrences.
Here, we need a Pair to indicate the string within the StringSet and a position within the string.

Write a program that traverses the nodes of the suffix tree of "mississippi" in the order shown here:

At each node print the text of the edges from the root to the node.
You may only use the functions goDown, goRight, goUp and isRoot to navigate and representative which returns the string that represents the node the iterator points to.

Modify the program to efficiently skip nodes with representatives longer than 3.
Move the whole program into a template function whose argument specifies the index type and call this function twice, once for the IndexEsa and once for the IndexWotd index.

Solution

We modify the DFS traversal to skip the descent if we walk into a node whose representative is longer than 3.
We then proceed to the right and up as long as the representative is longer than 3.

The FMIndex supports bidirectional iteration, i.e. a pattern can be extended to the left or right in an arbitrary order.
This is done by maintaining iterators on two separate indices, one on the original and one on the reversed text and keeping both iterators synchronized at all times.
The interface is similar to what you learned in the previous section.
All methods are extended by an additional tag specifying which iterator you want to use.
Going down the original iterator using the Fwd tag extends the pattern to the left (since the FMIndex is traversed as a prefix trie).
Using the Rev tag accesses the reversed text iterator and extends the pattern to the right.

Creating the index and iterator is very similar to unidirectional indices.
The FMIndex is wrapped in a BidirectionalIndex tag:

All methods for traversing the virtual trie are extended by the direction tag Fwd or Rev.
If none is used, it will access the iterator on the original text by default (same as using the Fwd tag).
The goUp method is the only method that does not specify a direction tag.
goUp corresponds to an undo operation, i.e. it rolls both iterators back to their previous states.

goDown(iter,DnaString("TTTC"),Fwd());// search CTTT in the prefix triegoDown(iter,Dna('G'),Rev());// extend to CTTTGgoUp(iter);std::cout<<representative(iter,Fwd())<<std::endl;std::cout<<representative(iter,Rev())<<std::endl;

Please bear in mind that you can also choose whether you want to retrieve the positions of hits in the original or reversed text:

// if we get here the pattern was found// output match positionsfor(unsignedi=0;i<length(getOccurrences(iter,Fwd()));++i)std::cout<<getOccurrences(iter,Fwd())[i]<<std::endl;for(unsignedi=0;i<length(getOccurrences(iter,Rev()));++i)std::cout<<getOccurrences(iter,Rev())[i]<<std::endl;

The tree traversal in assignment 2 is equal to a the tree traversal in a full depth-first search (dfs) over all suffix tree nodes beginning either in the root (preorder dfs) or in a leaf node (postorder dfs).
A preorder traversal (Preorder DFS) halts in a node when visiting it for the first time whereas a postorder traversal (Postorder DFS) halts when visiting a node for the last time.
The following two figures give an example in which order the tree nodes are visited.

Preorder DFS

Postorder DFS

Since these traversals are frequently needed SeqAn provides special iterators which we will describe next.

We want to construct the suffix tree of the string ‚Äúabracadabra‚ÄĚ and output the substrings represented by tree nodes in preorder dfs.
In order to do so, we create the string ‚Äúabracadabra‚ÄĚ and an index specialized with the type of this string.

#include<iostream>#include<seqan/index.h>usingnamespaceseqan;

The Iterator metafunction expects two arguments, the type of the container to be iterated and a specialization tag, as described earlier.
In this example we chose a TopDown History Iterator whose signature in the second template argument is TopDown<ParentLinks<Preorder>>.

If solely a postorder traversal is needed the BottomUp Iterator should be preferred as it is more memory efficient.
Please note that the BottomUp Iterator is only applicable to IndexEsa indices.

Tip

A relaxed suffix tree (see Indices) is a suffix tree after removing the $ characters and empty edges.
For some bottom-up algorithms it would be better not to remove empty edges and to have a one-to-one relationship between leaves and suffices.
In that cases you can use the tags PreorderEmptyEdges or PostorderEmptyEdges instead of Preorder or Postorder or EmptyEdges for the TopDown Iterator.

Note that the goNext is very handy as it simplifies the tree traversal in assignment 2 greatly.

Write a program that constructs an index of the StringSet ‚Äútobeornottobe‚ÄĚ, ‚Äúthebeeonthecomb‚ÄĚ, ‚Äúbeingjohnmalkovich‚ÄĚ and outputs the strings corresponding to suffix tree nodes in postorder DFS.

Solution

First we have to create a StringSet of CharString (shortcut for String<char>) and append the 3 strings to it.
This could also be done by using resize and then assigning the members with operator[].
The first template argument of the index class has to be adapted and is now a StringSet.

To switch to postorder DFS we have to change the specialization tag of ParentLinks from Preorder to Postorder.
Please note that the TopDownHistoryIterator always starts in the root node, which is the last postorder DFS node.
Therefore, the iterator has to be set explicitly to the first DFS node via goBegin.

Iterator<TMyIndex,TopDown<ParentLinks<Postorder>>>::TypemyIterator(myIndex);// Top-down iterators start in the root node which is not the first node of a// postorder DFS. Thus we have to manually go to the DFS start with goBegingoBegin(myIterator);while(!atEnd(myIterator)){std::cout<<representative(myIterator)<<std::endl;++myIterator;}

Alternatively with a TopDownHistoryIterator you also could have used a BottomUpIterator with the same result.
The BottomUp Iterator automatically starts in the first DFS node as it supports no random access.

As the last assignment lets try out one of the specialized iterators, which you can find at the bottom of this page.
Look there for the specialization which iterates over all maximal unique matches (MUMS).

After that we simply use the predefined iterator for searching MUMs, the MumsIterator.
Its constructor expects the index and optionally a minimum MUM length as a second parameter.
The set of all MUMs can be represented by a subset of suffix tree nodes.
The iterator will halt in every node that is a MUM of the minimum length.
The corresponding match is the node‚Äôs representative.

In the previous subsection we have seen how to walk through a suffix tree.
We now want to know what can be done with a suffix tree iterator.
As all iterators are specializations of the general VSTree Iterator class, they inherit all of its functions.
There are various functions to access the node the iterator points at (some we have already seen), so we concentrate on the most important ones.

returns the substring that represents the edge from the current node to its parent (only TopDownHistory Iterator)

Important

There is a difference between the functions isLeaf and isRightTerminal.
In a relaxed suffix tree (see Indices) a leaf is always a suffix, but not vice versa, as there can be internal nodes a suffix ends in.
For them isLeaf returns false and isRightTerminal returns true.

Some algorithms require to store auxiliary information (e.g. weights, scores) to the nodes of a suffix tree.
To attain this goal SeqAn provides so-called property maps, simple Strings of a property type.
Before storing a property value, these strings must first be resized with resizeVertexMap.
The property value can then be assigned or retrieved via assignProperty, getProperty, or property.
It is recommended to call resizeVertexMap prior to every call of assignProperty to ensure that the property map has sufficient size.
The following example iterates over all nodes in preorder dfs and recursively assigns the node depth to each node.
First we create a String of int to store the node depth for each suffix tree node.

The main loop iterates over all nodes in preorder DFS, i.e. parents are visited prior children.
The node depth for the root node is 0 and for all other nodes it is the parent node depth increased by 1.
The functions assignProperty, getProperty and property must be called with a VertexDescriptor.
The vertex descriptor of the iterator node is returned by value and the descriptor of the parent node is returned by nodeUp.

Given a string s a repeat is a substring r that occurs at 2 different positions i and j in s.
The repeat can also be identified by the triple (i,j,|r|).
A maximal repeat is a repeat that cannot be extended to the left or to the right, i.e. s[i-1]‚Č†s[j-1] and s[i+|r|]‚Č†s[j+|r|].
A supermaximal repeat r is a maximal repeat that is not part of another repeat.
Given a set of strings \(s_1, \dots, s_m\) a MultiMEM (multiple maximal exact match) is a substring r that occurs in
each sequence \(s_i\) at least once and cannot be extended to the left or to the right.
A MUM (maximal unique match) is a MultiMEM that occurs exactly once in each sequence.
The following examples demonstrate the usage of these iterators:

Indices in SeqAn allow efficient pattern queries in strings or sets of
strings. In contrast to, e.g., online-search algorithms that search
through the text in \(\mathcal{O}(n)\), substring indices find a
pattern in sublinear time \(o(n)\).

Indices store additional data structures that allow searching the text
using an iterator. Using the iterator can be thought of as traversing
a suffix tree. The following section gives you an introduction how
the suffix tree is built.

We consider an alphabet ő£ and a sentinel character $ that is smaller
than every character of ő£. A suffix tree of a given non-empty string s
over ő£ is a directed tree whose edges are labeled with non-empty
substrings of s$ with the following properties:

1. Each outgoing edge begins with a different letter and the outdegree
of an internal node is greater than 1.
2. Each suffix of s$ is the concatenation of edges from the root to a leaf node.
3. Each path from the root to a leaf node is a suffix of s$.

The following figure shows the suffix tree of the string s=‚ÄĚmississippi‚ÄĚ (suffix nodes are shaded):

Figure 1: Suffix tree of ‚Äúmississippi‚ÄĚ

Many suffix tree construction algorithms expect $ to be part of the
string alphabet which is undesirable for small bit-compressible
alphabets (e.g. DNA). In SeqAn there is no need to introduce a $. We
relax suffix tree criterion 2. and consider the relaxed suffix tree that
arises from the suffix tree of s by removing the $ character and all
empty edges. In the following, we only consider relaxed suffix trees and
simply call them suffix trees. In that tree a suffix can end in an inner
node as you can see in the next figure (suffix ‚Äúi‚ÄĚ):

This tutorial introduces you to the scoring systems that can be used in SeqAn to quantify the sequence similarity.
You will learn basic techniques to create and modify standard and custom scoring systems capable to satisfy the requirements of a wide range of applications.

The alignment procedures are usually based on the sequences similarity computation described by an alignment scoring system that gives countable information used to determine which sequences are related and which are not.

Four main biological events must be considered during the sequence alignment:
Conservation, substitution, insertion and deletion.
We could have a Conservation when the two compared letters are the same and a Match is detected, a Substitution when we detect a Mismatch where a letter is aligned with another, and Insertion or Deletion when in one of the two aligned sequences a letter is aligned with a Gap.
Matches, mismatches and gaps detected during the alignment do not guarantee to be the most representative biological truth since their dispositions is dependent of the chosen scoring schemes and the selected alignment algorithm. In order to improve the correlation between computed sequence alignment and biological similarity, specific combinations of scoring schemes and alignment algorithms have been developed during the years and are usually adopted for the alignment of different types of biological sequences. For example, as we will see in the following, the small RNA sequences are usually aligned with a Global Alignment algorithm implementing a Simple Score scheme, differently from the protein sequences that are mostly aligned with the Local Alignment algorithm that uses a Substitution Matrix Score scheme.

Given an alignment structure that store the two sequences and a scoring scheme, the score of the alignment can be computed as the sum of the scores for aligned character pairs plus the sum of the scores for all gaps.

With refer to the alignment procedure a Scoring Scheme can be defined as the set of rules used to assess the possible biological events that must be considered during the alignment procedure.

In SeqAn are available several scoring schemes to evaluate matches and mismatches, while three different gap models can be applied to consider insertions and deletions events.
We will first introduce you to the scoring schemes used to evaluate match and mismatch. Subsequently, you will learn how to chose the gap model to be implemented in the chosen scoring scheme.

The simplest example of Scoring Scheme, usually applied to score the similarity among nucleotide sequences, is the Levenshtein distance model that assigns a score of 0 and -1 respectively if a match or a mismatch occurs, whereas a penalty value equal to -1 in case of gaps representing insertions or deletions (this scoring scheme is the default for SimpleScore).
Alternatively, also the Hamming distance model can be used for some simple tasks that do not require the gap evaluations.

Now, let‚Äôs start by constructing our first scoring function for the global alignment algorithm called with the function globalAlignment.
As first step we need to include the header file <seqan/align.h> which contains the necessary data structures and functions associated with the alignments.
The next steps would be to implement the main function of our program and to define the types that we want to use.

We first define the type of the input sequences (TSequence) and an Align object (TAlign) type to store the alignment.
For more information on the Align datastructure, please read the tutorial Alignment Representation (Gaps).
After defining the types, we can continue to construct our own Align object.
First, we create two input sequences seq1="TELKDD" and seq2="LKTEL", then we define the scoring values for match, mismatch, gap.
As last we create the ‚Äėalign‚Äô object and resize it to manage two Gaps objects, at this point we filled it with the sequences to be aligned.

Now, we can compute the global alignment that makes use of the simple scoring function.
To do so, we simply call the function globalAlignment and give as input parameters the align object and the scoring scheme representing the Levenshtein distance.
The globalAlignment function fills the align object with the best computed alignment and returns the maximum score which we store in the score variable.
Afterwards, we print the computed score and the corresponding alignment.

Congratulations!
You have created your global alignment implementing the simple scoring function, the output is as follows:

Score: -5 0 . TELK-DD || --LKTEL

However, in the evaluation of protein similarity or for advanced nucleotide alignments a more complex scoring model is generally applied.
It is based on the usage of a Substitution Matrix, proven to better describe from a biological point of view, events such as matches and mismatches.

Substitutional Matrices are built on the basis of the probability that a particular amino acid or nucleotide is replaced with another during the evolution process.
They assign to each pair a value that indicates their degree of similarities, obtained thanks to statistical methods reflecting the frequency of a particular substitution in homologous protein or RNA families. A positive value in the Substitutional Matrix means that the two letters share identical or similar properties.

These scoring schemes store a score value for each pair of characters. This value can be accessed using score.
Examples for this kind of scoring scheme are Pam120 and Blosum62.
Anyway the class MatrixScore can be used to store arbitrary scoring matrices for the creation of custom scoring systems, as shown in the example proposed in the Working With Custom Score Matrices.

Blosum matrix, is one of the most used Substitutional Matrix implemented by considering multiple alignments of evolutionarily divergent proteins, while Ribosum is the RNA counterpart computed using ribosomal sequences.

In the following example it is proposed the construction of a scoring function for a global alignment algorithm that uses the Blosum62 matrix to score the matched and mismatched letters.
As first we include the header file <seqan/align.h> which contains the necessary data structures and functions associated with the alignments, then we implement the main function of our program and define the types that we want to use.

The input sequences type TSequence and the Align object of type TAlign are defined and the two input sequences seq1="TELKDD" and seq2="LKTEL" together with the gap penalty are assigned. In this case we define only the gap value since the Blosum matrix will be used to score matches and mismatches.
Then the sequences are associated with the alignment object.

Now, we compute the global alignment function, providing as second parameter the tag referred to the Blosum62 matrix together with the gap costs.
To do so, we simply call the function globalAlignment and give as input parameters the align object and the Blosum62 scoring scheme.
The globalAlignment function returns the score of the best alignment, which we store in the score variable that is then printed together with the corresponding alignment.

The output of a global alignment implementing the Blosum62 scoring function is as follows:

Score: 9 0 . --TELKDD ||| LKTEL---

Note

As can be noted the output of this scoring scheme is completely different with respect to the output generated with the simple scoring scheme confirming that the scoring scheme choice is one of the most important step to achieve high quality alignments.

In the previous sections we proposed two simple code examples useful to highlight the differences between two scoring schemes capable to evaluate match and mismatch events. In this section we will see the three gap models, implemented in the SeqAn library, to evaluate the insertion and deletion events.

The easiest is the Linear gap model that considers, for the alignment score computation, the gap length (g) giving the possibility to evaluate with different scores gaps of different sizes;

This gap model is chosen as standard when only a gap value is provided in the scoring function or when the two provided gaps have the same value. For instance, this gap model as been adopted during the alignment computation of the two proposed examples.

It has been proven that the first amino acid or nucleotide inserted/deleted (identified as gap open) found during the alignment operations is more significant, from a biological point of view, than the subsequent ones (called gap extension), making the so called Affine Gap model a viable solution for the alignment of biomolecules [Car06].
Affine gap model that attribute different costs to the gap open (d) and the gap extension (e) events, is able to assign an higher penalty to the gap presence with respect to its relative length (g).

The Affine Gap model implemented in the DP alignment algorithms is however quite expensive both in terms of computational time as well as in terms of memory requirements with respect to other less demanding solutions such as the Linear Gap model application.

In SeqAn is provided an optimised version of the Affine Gap model called Dynamic Gap Selector (DGS) designed by Urgese et al. [UPA+14]. This new gap model can be used to reduce the computational time and the memory requirement while keeping the alignment scores close to those computed with the Affine Gap model.
The usage of Dynamic Gap model in the Global alignment computation of long strings can give results slightly different from those computed using Affine Gap model since the alignment matrix became bigger and different alignment paths can be chosen during the alignment procedure. Score variation are rare when Dynamic Gap model is used in the Local alignments.

The order of the different costs in the scoring scheme is match, mismatch, gapExtend and gapOpen.
The gap model selection can be done providing one of the three specific tags (LinearGaps(), AffineGaps() or DynamicGaps()) as last parameter in the scoring function creation. If you want to use Linear Gap costs you could also omit the last parameter gapOpen and the scoring scheme would automatically choose the Linear Gap cost function.
The Affine Gap model is chosen as standard when the gap costs are different and the gap model tag is not provided. If the Dynamic Gap model is required the relative tag must be supplied.

In the following we propose an example where two different scoring functions have been created to show how to call a global alignment algorithm that uses the Blosum62 plus the AffineGaps() and DynamicGaps() specializations.
The inclusion of the header and the type definition is identical to the previous examples.

The input sequences type and the Align object of type TAlign are then create and initialized. As can be noted we define two different gap values, one for the gap extension and one for the gap open. Even in this example the Blosum62 will be used to score match and substitutions events.

Now, we can compute the global alignment function providing as second parameter the tag referred to the Blosum62 matrix filled with the two different gap costs. Moreover, the tag for the gap model selection is provided.
To do so, we simply call the function globalAlignment and give as input parameters the align object, the Blosum62 scoring scheme and the AffineGaps() or DynamicGaps() tag.
The globalAlignment function output is then printed.

This tutorial introduces you to the gaps data structures that can be used to represent an alignment in SeqAn.
You will learn basic techniques to create and modify such data structures and how to access certain information from these data structures.

The Align data structure is simply a set of multiple Gaps data structures.
A Gaps data structure is a container storing gap information for a given source sequence.
The gap information is put on top of the source sequence (coordinates of the gapped sequence refer to the gap space) without directly applying them to the source (coordinates of the ungapped sequence refer to the source space).
This way operating with gaps sustains very flexible.

There are two specializations available for the Gaps data structures:
Array Gaps and Anchor Gaps.
They differ in the way they implement the gap space.

Note

In general, using Array Gaps is sufficient for most applications.
This specialization is also the default one if nothing else is specified.
It simply uses an array which stores the counts of gaps and characters in an alternating order.
Thus, it is quite efficient to extend existing gaps while it is more expensive to search within the gapped sequence or insert new gaps.
Alternatively, one should prefer Anchor Gaps if many conversions between coordinates of the gap and the source space are needed as binary search can be conducted to search for specific positions.

Now, let‚Äôs start by constructing our first alignment.
Before we can make use of any of the mentioned data structures, we need to tell the program where to find the definitions.
This can be achieved by including the header file <seqan/align.h> which contains the necessary data structures and functions associated with the alignments.
The next steps would be to implement the main function of our program and to define the types that we want to use.

We first define the type of the input sequences (TSequence).
Then we can define the type of our actual Align object we want to use.
In an Align object, the gapped sequences are arranged in rows.
You can use the Metafunction Row to get the correct type of the used Gaps objects.
In the following we use the term row to explicitly refer to a single gapped sequence in the Align object.
We will use the term gappedsequence to describe functionalities that are related to the Gaps data structure independent of the Align object.

After defining the types, we can continue to actually construct our own Align object.
Therefore, we need to resize the alignment object in order to reserve space for the sequences we want to add.
In our case, we assume a pairwise alignment, hence we reserve space for 2 sequences.
With the function row, we get access to the gapped sequence at a specific row in the alignment object.
This is similar to the value function used in String Sets.
Now, we can assign the source to the corresponding gapped sequence.

The second string CDEFGAHGC of the alignment is cropped in the output to
CDEFGA, such that they are of equal length. Note that the string itself
is not modified, i.e. not shortened.

After assigning the sources to the gapped sequences, we need to add some gaps to make it look like a real alignment.
You can use the functions insertGap() and removeGap() to insert and delete one gap or insertGaps() and removeGaps() to insert and delete multiple gaps in a gapped sequence.

Congratulations!
You have created your first alignment.
Note that we used a reference declaration TRow& for the variables row1 and row2.
Without the reference, we would only modify copies of rows and the changes would not effect our align object.

In the next steps, we want to dig a little deeper to get a feeling for the gap space and the source space.
As mentioned above, the gaps are not inserted into the source but put on top of it in a separate space, the gap space.
When inserting gaps, the gap space is modified and all coordinates right of the inserted gap are shifted to the right by the size of the gap.
At the same time, the coordinates of the source remain unchanged.
Using the function toSourcePosition(), we can determine which position of the current gapped sequence (gap space) corresponds to the position in the source space.

If the position in the gap space is actually a gap, then toSourcePosition() returns the source position of the next character to the right that is not a gap.
Vice versa, we can determine where our current source position maps into the gap space using the function toViewPosition().

In the first alignment, it seems that the end of the second row is cropped off to match the size of the first one.
This effect takes place only in the visualization but is not explicitly applied to the gapped sequence.
The second alignment is the one we manually constructed.
Here, you can see that the second row is expanded to its full size while it matches the size of the first row.
However, it is possible to explicitly crop off the ends of a gapped sequence by using the functions setClippedBeginPosition() and setClippedEndPosition().
These functions shrink the gap space and can be understood as defining an infix of the gapped sequence.
After the clipping, the relative view position changes according to the clipping and so does the mapping of the source positions to the gap space.
The mapping of the view positions to the source space does not change.

It is important to understand the nature of the clipping information.
It virtually shrinks the gap space not physically.
That means the information before/after the begin/end of the clipping still exists and the physical gap space remains unchanged.
To the outer world it seems the alignment is cropped off irreparably.
But you can expand the alignment again by resetting the clipping information.

In the last part of this section, we are going to iterate over a Gaps object.
This can be quite useful if one needs to parse the alignment rows to access position specific information.
First, we have to define the type of the Iterator, which can be easily done by using the metafunction Iterator.
Remember that we iterate over an TRow object.
Then we have to construct the iterators it which points to the begin of row1 using the begin() function and itEnd which points behind the last value of row1 using the end() function.
If you need to refresh the Iterator Concept you can read the iterator section Iteration.
While we iterate over the gapped sequence, we can ask if the current value, at which the iterator it points to, is a gap or not by using the function isGap().
Use gapValue to print the correct gap symbol.

Construct an alignment using the Align data structure for the sequences "ACGTCACCTC" and "ACGGGCCTATC".
Insert two gaps at the second position and insert one gap at the fifth position of the first sequence.
Insert one gap at the ninth position of the second sequence.
Iterate over the rows of your Align object and print the total count of gaps that exist within the alignment.

Hints

You can use the function countGaps to count the number of consecutive gaps starting from the current position of the iterator.

The resulting alignment should look like:

AC--GTC-ACCTCACGGGCCTA--TC

Solution

#include<iostream>#include<seqan/align.h>usingnamespaceseqan;intmain(){// Defining all types that are needed.typedefString<char>TSequence;typedefAlign<TSequence,ArrayGaps>TAlign;typedefRow<TAlign>::TypeTRow;typedefIterator<TRow>::TypeTRowIterator;TSequenceseq1="ACGTCACCTC";TSequenceseq2="ACGGGCCTATC";// Initializing the align object.TAlignalign;resize(rows(align),2);assignSource(row(align,0),seq1);assignSource(row(align,1),seq2);// Use references to the rows of align.TRow&row1=row(align,0);TRow&row2=row(align,1);// Insert gaps.insertGaps(row1,2,2);insertGap(row1,7);// We need to pass the view position which is changed due to the previous insertion.insertGaps(row2,9,2);// Initialize the row iterators.TRowIteratoritRow1=begin(row1);TRowIteratoritEndRow1=end(row1);TRowIteratoritRow2=begin(row2);// Iterate over both rows simultaneously.intgapCount=0;for(;itRow1!=itEndRow1;++itRow1,++itRow2){if(isGap(itRow1)){gapCount+=countGaps(itRow1);// Count the number of consecutive gaps from the current position in row1.itRow1+=countGaps(itRow1);// Jump to next position to check for gaps.itRow2+=countGaps(itRow1);// Jump to next position to check for gaps.}if(isGap(itRow2)){gapCount+=countGaps(itRow2);// Count the number of consecutive gaps from the current position in row2.itRow1+=countGaps(itRow2);// Jump to next position to check for gaps.itRow2+=countGaps(itRow2);// Jump to next position to check for gaps.}}// Print the result.std::cout<<"Number of gaps: "<<gapCount<<std::endl;}

This tutorial introduces you to the graph data structures that can be used to represent an alignment in SeqAn.
You will learn basic techniques to create and modify such data structures and how to access certain information from these data structures.

Another very useful representation of alignments is given by the Alignment Graph.
It is a graph in which each vertex corresponds to a sequence segment, and each edge indicates an ungapped alignment between the connected vertices, or more precisely between the sequences stored in those vertices.
Here is an example of such a graph:

In the following we will actually construct this example step by step.
First we include the iostream header from the STL and the <seqan/align.h> header to include all necessary functions and data structures we want to use.
We use the namespace seqan and write the main function with an empty body.

At the begin of the function we define our types we want to use later on.
We define TSequence as the type of our input strings.
Since we work with a Dna alphabet we define TSequence as a String over a Dna alphabet.
For the AlignmentGraph we need two StringSets.
The TStringSet is used to actually store the input sequences and the TDepStringSet is internally used by the AlignmentGraph.
That is the AlignmentGraph does not copy the sources into its data structure but rather stores a reference to each of the given input strings as it does not modify the input sequences.
The Dependent StringSet facilitates this behavior.
In the end we define the actual AlignmentGraph type.

We first create our two input sequences TTGT and TTAGT append them to the StringSet strings using the appendValue function and pass the initialized strings object as a parameter to the constructor of the AlignmentGraph alignG.

Before adding vertices to the graph align prints the empty adjacency and edge
list.

Adjacency list:Edge list:

Before we construct the alignment we print the unmodified AlignmentGraph.
Then we add some alignment information to the graph.
In order to add an ungapped alignment segment we have to add an edge connecting two vertices of different input sequences.
To do so we can use the function addEdge and specify the two vertices that should be connected.
Since we do not have any vertices yet, we create them on the fly using the function addVertex().
The function addVertex gets as second parameter the id which points to the the correct input sequence within the strings object.
We can use the function positionToId() to receive the id that corresponds to a certain position within the underlying Dependent StringSet of the AlignmentGraph.

We can access the Dependent StringSet using the function stringSet().
The third parameter of addVertex specifies the begin position of the segment within the respective input sequence and the fourth parameter specifies its length.
Now, we add an edge between the two vertices of each input sequence which covers the first two positions.
In the next step we have to add a gap.
We can do this simply by just adding a vertex that covers the inserted string.
Finally we have to add the second edge to represent the last ungapped sequence and then we print the constructed alignment.

Note that we use findVertex() to find the the
last two inserted vertices. The syntax is the same as addVertex(), but omits the length parameter.

Construct a multiple sequence alignment using the Alignment Graph data structure.
Use the three sequences GARFIELDTHECAT, GARFIELDTHEBIGCAT and THEBIGCAT and align them such that you obtain the maximal number of matches.

Hints

TSequence should be String<char> instead of String<Dna>.

The function findVertex returns the vertex of an AlignmentGraph that covers the given position in the given sequence.

Solution

#include<iostream>#include<seqan/align.h>usingnamespaceseqan;intmain(){// Define the types we need.typedefString<char>TSequence;typedefStringSet<TSequence>TStringSet;typedefStringSet<TSequence,Dependent<>>TDepStringSet;typedefGraph<Alignment<TDepStringSet>>TAlignGraph;typedeftypenameVertexDescriptor<TAlignGraph>::TypeTVertexDescriptor;// Initializing the sequences and the string set.TSequenceseq1="GARFIELDTHECAT";TSequenceseq2="GARFIELDTHEBIGCAT";TSequenceseq3="THEBIGCAT";TStringSetstrings;appendValue(strings,seq1);appendValue(strings,seq2);appendValue(strings,seq3);// Load the string set into the Alignment Graph.TAlignGraphalignG(strings);TVertexDescriptoru,v;// Add two vertices covering "GARFIELD" in the first and the second sequence and connect them with an edge.u=addVertex(alignG,positionToId(stringSet(alignG),0),0,8);v=addVertex(alignG,positionToId(stringSet(alignG),1),0,8);addEdge(alignG,u,v);// Add two vertices covering "THE" in the first and the second sequence and connect them with an edge.u=addVertex(alignG,positionToId(stringSet(alignG),0),8,3);v=addVertex(alignG,positionToId(stringSet(alignG),1),8,3);addEdge(alignG,u,v);// Find the vertex covering "THE" in the first sequence and add the vertex covering "THE" in the third sequence and connect them with an edge.u=findVertex(alignG,positionToId(stringSet(alignG),0),8);v=addVertex(alignG,positionToId(stringSet(alignG),2),0,3);addEdge(alignG,u,v);// Find the vertices covering "THE" in the second and the third sequence and connect them with an edge.u=findVertex(alignG,positionToId(stringSet(alignG),1),8);v=findVertex(alignG,positionToId(stringSet(alignG),2),0);addEdge(alignG,u,v);// Add two vertices covering "FAT" in the second and the third sequence and connect them with an edge.u=addVertex(alignG,positionToId(stringSet(alignG),1),11,3);v=addVertex(alignG,positionToId(stringSet(alignG),2),3,3);addEdge(alignG,u,v);// Add two vertices covering "CAT" in the first and the second sequence and connect them with an edge.u=addVertex(alignG,positionToId(stringSet(alignG),0),11,3);v=addVertex(alignG,positionToId(stringSet(alignG),1),14,3);addEdge(alignG,u,v);// Find the vertex covering "CAT" in the first sequence and add the vertex covering "CAT" in the third sequence and connect them with an edge.u=findVertex(alignG,positionToId(stringSet(alignG),0),11);v=addVertex(alignG,positionToId(stringSet(alignG),2),6,3);addEdge(alignG,u,v);// Find the vertices covering "CAT" in the second and the third sequence and connect them with an edge.u=findVertex(alignG,positionToId(stringSet(alignG),1),14);v=findVertex(alignG,positionToId(stringSet(alignG),2),6);addEdge(alignG,u,v);std::cout<<alignG<<std::endl;return0;}

Alignment Algorithms ( e.g.
Pairwise and
Multiple )
are one of the core algorithms in SeqAn. In this section you can learn how SeqAn
represents alignments as C++ Objects and how you could use those data structures
for your own alignment algorithm. Furthermore, you can learn which different
kinds of Scoring Schemes exist,
i.e. which combinations of Match/Mismatch Evaluation (e.g. Simple Score,
Substitutional Matrices Score) and Insertion/Deletion Evaluation (e.g. Linear
Gap Model, Affine Gap Model and Dynamic Gap Model) are possible and how you can
define your own Scoring Matrices.

You will learn how to work with annotations in SeqAn.
After this tutorial, you will be able to write your own programs using annotations and analyzing them.
You will be ready to continue with the Fragment Store Tutorial, e.g. if you want to combine your annotations with information from alignments.

This tutorial will present SeqAn‚Äôs efficient and easy-to-use data structures to work with annotations.
They allow to annotate genome regions with features like ‚Äėgene‚Äô, ‚ÄėmRNA‚Äô, ‚Äėexon‚Äô, ‚Äėintron‚Äô and if required with custom features.
We will give you an understanding of how to load annotations from a GFF or GTF file, store them in efficient data structures, as well as how to traverse and access these information.

This section will give you a short introduction to data structures relevant for working with annotations.

In SeqAn, annotations are stored in the so-called annotationStore, which is part of the FragmentStore.
The annotationStore can only be used together with the FragmentStore, because the latter stores additional information, e.g. the contig names or sequences.
The FragmentStore is a data structure specifically designed for read mapping, genome assembly or gene annotation.

The FragmentStore can be seen as a database, where each table (called ‚Äústore‚ÄĚ) is implemented as a String.
Each row of the table corresponds to an element in the string.
The position of each element in the string implicitly represents the Id of such element in the table.
All such strings are members of the class FragmentStore, are always present and empty if unused.
For example, the member contigStore is a string of elements, each one containing among others a contig sequence.

Before we deal with the actual annotation tree, we will first describe how you can easily load annotations from a GFF or GTF file into the FragmentStore.

An annotation file can be read from an GffFileIn with the function readRecords.
The file extension specifies if we want to read a GFF, GTF or UCSC file.
The following example shows how to read an GTF file:

The GFF-reader is also able to detect and read GTF files.
The UCSC Genome Browser uses two separate files, the knownGene.txt and knownIsoforms.txt.
They must be read by using two different UcscFileIn objects (one for knownGene.txt and one for knownIsoforms.txt).
Finally you call readRecords with both UcscFileIn objects.

Tip

An annotation can be loaded without loading the corresponding contigs.

In that case empty contigs are created in the contigStore with names given in the annotation.
A subsequent call of loadContigs would load the sequences of these contigs, if they have the same identifier in the contig file.

This section will illustrate how to use iterators to traverse the annotation tree.

The annotation tree can be traversed and accessed with the AnnotationTree Iterator.
Again we use the metafunction Iterator to determine the appropriate iterator type for our container.
A new AnnotationTree iterator can be obtained by calling begin with a reference to the FragmentStore and the AnnotationTree tag:

The AnnotationTree iterator starts at the root node and can be moved to adjacent tree nodes with the functions goDown, goUp, and goRight.
These functions return a boolean value that indicates whether the iterator could be moved.
The functions isLeaf, isRoot, isLastChild return the same boolean without moving the iterator.
With goRoot or goTo the iterator can be moved to the root node or an arbitrary node given its annotationId.
If the iterator should not be moved but a new iterator at an adjacent node is required, the functions nodeDown, nodeUp, nodeRight can be used.

// Move the iterator down to a leafwhile(goDown(it)){}// Create a new iterator and if possible move it to the right sibling of the first iteratorIterator<FragmentStore<>,AnnotationTree<>>::Typeit2;if(isLastChild(it))it2=nodeRight(it);

The AnnotationTree iterator supports a preorder DFS traversal and therefore can also be used in typical begin-end loops with the functions goBegin (== goRoot), goEnd, goNext, atBegin, atEnd.
During a preorder DFS, the descent into subtree can be skipped by goNextRight, or goNextUp which proceeds with the next sibling or returns to the parent node and proceeds with the next node in preorder DFS.

// Move the iterator back to the beginninggoBegin(it);// Iterate over the nodes in preorder DFS while the end is not reached and// output if the current node is a leafwhile(!atEnd(it)){if(isLeaf(it))std::cout<<" current node is leaf"<<std::endl;goNext(it);}

Copy the code below, which loads the annotations from a given GTF file into the FragmentStore and initializes an iterator on the AnnotationTree.
Download the GTF file assignment_annotations.gtf, whose annotations build an AnnotationTree of the typical structure with gene, mRNA and exon level.
Adjust the code to go down to the exon level and iteratate over all children of the first mRNA and count them.
Print the result.

In the given data the left-most leaf is a child of mRNA and has siblings.
You can use the function goRight to traverse over all siblings.

Solution

Click more‚Ä¶ to see one possible solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>#include<seqan/store.h>usingnamespaceseqan;intmain(){CharStringfileName=getAbsolutePath("demos/tutorial/genome_annotations/assignment_annotations.gtf");GffFileInfile(toCString(fileName));FragmentStore<>store;readRecords(store,file);// Create AnnotationTree iteratorIterator<FragmentStore<>,AnnotationTree<>>::Typeit;it=begin(store,AnnotationTree<>());unsignedcount=0;// Go down to the first leaf (first child of the first mRNA)while(goDown(it)){}std::cout<<"Is leaf: "<<isLeaf(it)<<std::endl;++count;// Iterate over all siblings and countwhile(goRight(it))++count;std::cout<<"No. of children of the first mRNA: "<<count<<std::endl;return0;}

Reuse the code and the GTF file from above.
Instead of counting only the children of the first mRNA adjust the code to count the children for each given mRNA.
Print the results.

Hints

After you reached the last child of the first mRNA you can use the functions goNext and goDown to traverse to the next leaf.

Solution

Click more‚Ä¶ to see one possible solution.

#include<iostream>#include<seqan/sequence.h>#include<seqan/stream.h>#include<seqan/store.h>usingnamespaceseqan;intmain(){CharStringfileName=getAbsolutePath("demos/tutorial/genome_annotations/assignment_annotations.gtf");GffFileInfile(toCString(fileName));FragmentStore<>store;readRecords(store,file);// Iterate over all leafs, count and print the resultIterator<FragmentStore<>,AnnotationTree<>>::Typeit;it=begin(store,AnnotationTree<>());unsignedcount=0;std::cout<<"Number of children for each mRNA: "<<std::endl;// Go down to the first leaf (first child of the first mRNA)while(goDown(it)){}while(!atEnd(it)){++count;// Iterate over all siblings and countwhile(goRight(it))++count;std::cout<<count<<std::endl;count=0;// Jump to the next mRNA or gene, go down to its first leaf and count itif(!atEnd(it)){goNext(it);if(!atEnd(it))while(goDown(it)){}}}return0;}

Let us now have a closer look how to access the information stored in the different stores representing the annotation tree.

To access or modify the node an iterator points at, the iterator returns the node‚Äôs annotationId by the value function (== operator*).
With the annotationId the corresponding entry in the annotationStore could be modified manually or by using convenience functions.
The function getAnnotation returns a reference to the corresponding entry in the annotationStore.
getName and setName can be used to retrieve or change the identifier of the annotation element.
As some annotation file formats don‚Äôt give every annotation a name, the function getUniqueName returns the name if non-empty or generates one using the type and id.
The name of the parent node in the tree can be determined with getParentName.
The name of the annotation type, e.g. ‚ÄėmRNA‚Äô or ‚Äėexon‚Äô, can be determined and modified with getType and setType.

Assume we have loaded the file example.gtf with the following content to the FragmentStorestore and instantiated the iterator it of the corresponding annotation tree.

We now want to iterate to the first exon and output a few pieces of information:

// Move the iterator to the begin of the annotation treeit=begin(store,AnnotationTree<>());// Go down to exon levelwhile(goDown(it));std::cout<<"type: "<<getType(it)<<std::endl;std::cout<<"id: "<<value(it)<<std::endl;std::cout<<"begin position: "<<getAnnotation(it).beginPos<<std::endl;

For our example the output would be:

type: exonid: 3begin position: 149

An annotation can not only refer to a region of a contig but also contain additional information given as key-value pairs.
The value of a key can be retrieved or set by getValueByKey and assignValueByKey.
The values of a node can be cleared with clearValues.

You will learn about the SeqAn FragmentStore for handling fragments.
The ‚Äúfragments‚ÄĚ are reads and the data structure is useful in the context of read mapping, genome assembly, and gene annotation.
After completing this tutorial, you will be able to use the most relevant functionality of the FragmentStore class.

The FragmentStore is a data structure specifically designed for read mapping, genome assembly or gene annotation.
These tasks typically require lots of data structures that are related to each other like:

reads, mate-pairs, reference genome

pairwise alignments

genome annotation

The Fragment Store subsumes all these data structures in an easy to use interface.
It represents a multiple alignment of millions of reads or mate-pairs against a reference genome consisting of multiple contigs.
Additionally, regions of the reference genome can be annotated with features like ‚Äėgene‚Äô, ‚ÄėmRNA‚Äô, ‚Äėexon‚Äô, ‚Äėintron‚Äô or custom features.
The Fragment Store supports I/O functions to read/write a read alignment in SAM/BAM or AMOS format and to read/write annotations in GFF or GTF format.

The Fragment Store can be compared with a database where each table (called ‚Äústore‚ÄĚ) is implemented as a String member of the FragmentStore class.
The rows of each table (implemented as structs) are referred by their ids which are their positions in the string and not stored explicitly (marked with * in the Figures 2 and 5).
The only exception is the alignedReadStore whose elements of type AlignedReadStoreElement contain an id-member as they may be rearranged in arbitrary order, e.g. by increasing genomic positions or by readId.
Many stores have an associated name store to store element names.
Each name store is a StringSet that stores the element name at the position of its id.
All stores are present in the Fragment Store and empty if unused.
The concrete types, e.g. the position types or read/contig alphabet, can be easily changed by defining a custom config struct which is a template parameter of the Fragment Store class.

The Fragment Store can represent a multiple read alignment, i.e. is an alignment between the contigs and the set of reads, where one read can be aligned at zero, one or multiple positions of a contig.
In the multiple alignment the contig is represented by one line with gaps (-) and the remaining lines are to reads or read segments with gaps aligned to the contig.
The following figure shows one contig (the line at the top) and multiple reads aligned to it arranged as stairs (reads in lower-case align to the reverse strand):

The following figure shows which tables represent the multiple read alignment:

*Figure 2:* Stores used to represent a multiple read alignment

The main table is the alignedReadStore which stores AlignedReadStoreElements.
Each entry is an alignment of a read (readId) and a contig (contigId).
Introduced gaps are stored as a string of gap anchors in the gaps member of the alignedReadStore entry and the contigStore entry.
The begin and end positions of the alignment are given by the beginPos and endPos members which are 0-based positions on the forward strand in gap space, i.e. positions in the gapped contig sequence.
If the read is aligned to the reverse strand it holds endPos<beginPos.
However, the gaps are always related to the forward strand.
Additional information, e.g. the number of errors, an alignment score or additional alignment tags, are stored in the tables alignQualityStore and alignedReadTagStore at position id, where id is a unique id of the AlignedReadStoreElement.
Paired-end or mate pair alignments are represented by two entries in the alignedReadStore that have the same pairMatchId value (unequal to INVALID_ID).
For orphaned read alignments holds pairMatchId==INVALID_ID.

The alignedReadStore is the only store where the id (alignId in the figure) of an element is not implicitly given by its position.
The reason for this is that it is necessary in many cases to rearrange the elements of the alignedReadStore, e.g. increasingly by (contigId,beginPos), by readId or pairMatchId.
This can be done by sortAlignedReads.
If it is necessary to address an element by its id, the elements must be sorted by id first.
In the case that ids are not contiguously increasing, e.g. because some elements where removed, they must be renamed by a prior call of compactAlignedReads.
Analogously the function compactPairMatchIds renames pairMatchId values contiguously and replaces values that occur in only one alignment by INVALID_ID.

The multiple read alignment can be displayed in text form or in a scalable graphics format (SVG).
Therefore first a stairs layout of the reads must be computed via layoutAlignment and stored in an AlignedReadLayout.
The function printAlignment can then be used to output a window (beginPos,endPos,firstLine,lastLine) of the read alignment against a contig either to a stream or SVGFile.
The following small example demonstrates how to first load two contigs from a Fasta file and then import a read alignment given in SAM format:

In the next step, we want to access several pairwise alignments between reads and contig segments.
Therefore we first need to get the associated types that the Fragment Store uses to store contig and read sequences and gaps.
This can be done by the following typedefs:

Now we want to extract and output the alignments from the alignedReadStore at position 140,144,‚Ä¶,156.
First we store a reference of the alignedRead in ar as we need to access it multiple times.
The read sequence is neither stored in the readStore or alignedReadStore as many short sequences can more efficiently be stored in a separate StringSet like the readSeqStore.
We copy the read sequence into a local variable (defined outside the loop to save allocations/deallocations) as we need to compute the reverse-complement for reads that align to the reverse strand.
Then we create a gaps object that represent the alignment rows of the contig and the aligned read in the multiple sequence alignment.
The Gaps object requires references of the sequence and the gap-anchor string stored in the contigStore and the alignedReadStore.
We need to limit the view of the contig alignment row to the interval the read aligns to, i.e. the gap position interval [beginPos,endPos[.
After that we output both alignment rows.

Tip

The Gaps contains two Holder references to the sequence and the inserted gaps.
In our example these Holders are dependent and changes made to the Gaps object like the insertion/deletion of gaps would immediatly be persistent in the Fragment Store.

Then, we not only need to reverse-complement readSeq if the read aligns to the reverse strand (endPos<beginPos) but also need to convert its letters into lower-case.
Therefor SeqAn provides the function toLower.
Alternatively, we could iterate over readSeq and add (‚Äėa‚Äô-‚ÄėA‚Äô) to its elements.

It starts at the root node and can be moved to adjacent tree nodes with the functions goDown, goUp, and goRight.
These functions return a boolean value that indicates whether the iterator could be moved.
The functions isLeaf, isRoot, isLastChild return the same boolean without moving the iterator.
With goRoot or goTo it can be moved to the root node or an arbitrary node given its annotationId.
If the iterator should not be moved but a new iterator at an adjacent nodes is required, the functions nodeDown, nodeUp, nodeRight can be used.

The AnnotationTree iterator supports a preorder DFS traversal and therefore can also be used in typical begin-end loops with the functions goBegin (== goRoot), goEnd, goNext, atBegin, atEnd.
During a preorder DFS, the descent into subtree can be skipped by goNextRight, or goNextUp which proceeds with the next sibling or returns to the parent node and proceeds with the next node in preorder DFS.

To access or modify the node an iterator points at, the iterator returns the node‚Äôs annotationId by the value function (== operator*).
With the annotationId the corresponding entry in the annotationStore could be modified manually or by using convenience functions.
The function getAnnotation returns a reference to the corresponding entry in the annotationStore.
getName and setName can be used to retrieve or change the identifier of the annotation element.
As some annotation file formats don‚Äôt give every annotation a name, the function getUniqueName returns the name if non-empty or generates one using the type and id. The name of the parent node in the tree can be determined with getParentName.
The name of the annotation type, e.g. ‚ÄėmRNA‚Äô or ‚Äėexon‚Äô, can be determined and modified with getType and setType.

An annotation can not only reference a region of a contig but also contain additional information given as key-value pairs.
The value of a key can be retrieved or set by getValueByKey and assignValueByKeq.
The values of a node can be cleared with clearValues.

To efficiently load reads, use the function loadReads which auto-detects the file format, supporting Fasta, Fastq, QSeq and Raw (see AutoSeqFormat), and uses memory mapping to efficiently load millions of reads, their names and quality values.
If not only one but two file names are given, loadReads loads mate pairs or paired-end reads stored in two separate files.
Both files are required to contain the same number or reads and reads stored at the same line in both files are interpreted as pairs.
The function internally uses appendRead or appendMatePair and reads distributed over multiple files can be loaded with consecutive calls of loadReads.

Contigs can be loaded with the function loadContigs.
The function loads all contigs given in a single file or multiple files given a single file name or a StringSet of file names.
The function has an additional boolean parameter loadSeqs to load immediately load the contig sequence or if false load the sequence later with loadContig to save memory, given the corresponding contigId.
If the contig is accessed by multiple instances/threads the functions lockContig and unlockContig can be used to ensure that the contig is loaded and release it after use.
The function unlockAndFreeContig can be used to clear the contig sequence and save memory if the contig is not locked by any instance.

As SAM supports a multiple read alignment (with padding operations in the CIGAR string) but does not enforce its use.
That means that a typical SAM file represents a set of pairwise (not multiple) alignments.
To convert all the pairwise alignments into a multiple alignments of all reads, read internally calls the function convertPairWiseToGlobalAlignment.
A prior call to loadReads is not necessary (but possible) as SAM contains the read names, sequences and quality values.
Contigs can be loaded at any time.
If they are not loaded before reading a SAM file, empty sequences are created with the names referred in the SAM file.
A subsequent call of loadContigs would load the sequences of these contigs, if they have the same identifier in the contig file.

The GffFileIn is also able to detect and read GTF files in addition to GFF files.
As the knownGene.txt and knownIsoforms.txt files are two separate files used by the UCSC Genome Browser, they must be read by two consecutive calls of readRecords (first knownGene.txt then knownIsoforms.txt).
An annotation can be loaded without loading the corresponding contigs.
In that case empty contigs are created in the contigStore with names given in the annonation.
A subsequent call of loadContigs would load the sequences of these contigs, if they have the same identifier in the contig file.

Please note, that UCSC files cannot be written due to limitations of the file format.

In this tutorial we will explain the GenomeAnnotation store and the FragmentStore.
Both are very comprehensive and feature rich data structures that come with many advantages when dealing with annotation trees or read alignments.
However, these are more advanced topics and rather aimed for more experienced SeqAn users.

A graph in computer science is an ordered pair \(G = (V, E)\) of a set of vertices V and a set of edges E.
SeqAn provides different types of graphs and the most well-known graph algorithms as well as some specialized alignment graph algorithms.
In this part of the tutorial, we demonstrate how to construct a graph in SeqAn and show the usage of some algorithms.
Alignment graphs are described in the tutorial Alignment.

Let us follow a simple example.
We have given the following network of five cities and in this network we want to compute the shortest path from Hannover to any other city.

In the section Graph Basics, we will create the network and write the graph to a .dot file.
The section Property Maps assigns city names to the vertices and Graph Iterators demonstrates the usage of a vertex iterator.

After having worked through these sections you should be familiar with the general usage of graphs in SeqAn. You are then prepared to proceed with Graph Algorithms, where we will compute the shortest path by calling a single function.

The general header file for all types of graphs is <seqan/graph_types.h>.
It comprises the Graph class, its specializations, every function for basic graph operations and different iterators.
Later, for computing the shortest path we will also need <seqan/graph_algorithms.h> which includes the implementations of most of SeqAn‚Äôs graph algorithms.

We want to model the network of cities as an undirected graph and label the edges with distances.
Before we start creating edges and vertices, we need some typedefs to specify the graph type.

SeqAn offers different specializations of the class Graph: Undirected Graph, DirectedGraph, Tree, Word Graph, Automaton, HmmGraph, and Alignment Graph.
For our example, an undirected graph will be sufficient, so we define our own graph type TGraph with the specialization Undirected Graph of the class Graph.
Luckily, this specialization has an optional cargo template argument, which attaches any kind of object to the edges of the graph.
This enables us to store the distances between the cities, our edge labels, using the cargo type TCargo defined as unsignedint.
Using the cargo argument, we have to provide a distance when adding an edge.
And when we remove an edge we also remove the distance.

Each vertex and each edge in a graph is identified by a so-called descriptor.
The type of the descriptors is returned by the metafunction VertexDescriptor.
In our example, we define a type TVertexDescriptor by calling VertexDescriptor on our graph type.
Analogously, there is the metafunction EdgeDescriptor for edge descriptors.

We can now create the graph g of our type TGraph.

TGraphg;

For our example, we add five vertices for the five cities, and six edges connecting the cities.

Vertices can be added to g by a call to the function addVertex.
The function returns the descriptor of the created vertex.
These descriptors are needed to add the edges afterwards.

The function addEdge adds an edge to the graph.
The arguments of this function are the graph to which the edge is added, the vertices that it connects, and the cargo (which is in our case the distance between the two cities).

The function addEdges takes as parameters an array of vertex descriptors and the number of edges.
The array of vertex descriptors is sorted in the way predecessor1, successor1, predecessor2, successor2, ‚Ä¶