SeqAn Manual¶

Library Design Aims¶

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.

1. 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.
2. Consistency. Be consistent wherever possible, even at slight costs of efficiency.
3. Ease of use. The library should be easily usable wherever possible, even at slight costs of efficiency.
4. Reusability and Generosity. The algorithms in SeqAn should be reusable and generic, even at small costs of efficiency.

Modern C++¶

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++).

References

References are alternatives to pointers in C++ to construct value aliases. Also see Wikipedia on C++ references.

Templates

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 ? extends T 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.

Memory Management in SeqAn¶

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.

#include <seqan/sequence.h>

int main(int argc, char const ** argv)
{
    seqan::String<char> programName = argv[0];
    if (argc > 1)
    {
        seqan::String<char> firstArg = argv[1];
        if (argc > 2)
            return 1;
    }
    return 0;
}

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 return 0. 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 return 0 or the return 1. Corresponding code in C would be (arguably) more messy, either requiring goto or multiple free() calls, one before either return.

Motivation for Template Programming¶

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.

OOP vs. Generic Programming¶

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 Steps Tutorial.

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 :)!

Global Function Interface¶

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 get in touch with 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 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!

Meta-Programming¶

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 task 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.

    String<AminoAcid> amino_str = "ARN";

void amino_exchangeFirstValues(String<AminoAcid> & str)
{
    if (length(str) < 2)
        return;
    AminoAcid temp = str[0];
    str[0] = str[1];
    str[1] = temp;
}

template <typename T>
void general_exchangeFirstValues(T & str)
{
    if (length(str) < 2)
        return;
    AminoAcid temp = str[0];
    str[0] = str[1];
    str[1] = temp;
}

template <typename T>
void exchangeFirstValues(T & str)
{
    if (length(str) < 2)
        return;
    typename Value<T>::Type temp = str[0];
    str[0] = str[1];
    str[1] = temp;
}

And now?¶

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!

Running Example¶

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 following figure shows an expected result:

score:    101 ...        ... 801 ...
text:     This is an awesome tutorial to get to know SeqAn!
pattern:  tutorial           tutorial
           tutorial           tutorial
            ...                ...

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>

using namespace seqan;

int main()
{
    // Initialization
    String<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 (unsigned i = 0; i < length(text) - length(pattern) + 1; ++i)
    {
        int localScore = 0;
        // Iteration over the pattern for character comparison
        for (unsigned j = 0; j < length(pattern); ++j)
        {
            if (text[i + j] == pattern[j])
                ++localScore;
        }
        score[i] = localScore;
    }

    // Printing the result
    for (unsigned i = 0; i < length(score); ++i)
        std::cout << score[i] << " ";
    std::cout << std::endl;

    return 0;
}

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 using namespace seqan; 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.

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

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.

SeqAn and Templates¶

Let us now have a detailed look at the program.

We first include the IOStreams library that we need to print to the screen and the SeqAn’s <seqan/file.h> as well as <seqan/sequence.h> module from the SeqAn library that provides SeqAn String.

#include <iostream>
#include <seqan/file.h>
#include <seqan/sequence.h>

using namespace seqan;

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.

    // Initialization
    String<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);

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 (unsigned i = 0; i < length(text) - length(pattern) + 1; ++i)
    {
        int localScore = 0;
        // Iteration over the pattern for character comparison
        for (unsigned j = 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.

    // Printing the result
    for (unsigned i = 0; i < length(score); ++i)
        std::cout << score[i] << " ";
    std::cout << std::endl;

Refactoring¶

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>

using namespace seqan;

int computeLocalScore(String<char> subText, String<char> pattern)
{
    int localScore = 0;
    for (unsigned i = 0; i < length(pattern); ++i)
        if (subText[i] == pattern[i])
            ++localScore;

    return localScore;
}

String<int> computeScore(String<char> text, String<char> pattern)
{
    String<int> score;
    resize(score, length(text) - length(pattern) + 1, 0);

    for (unsigned i = 0; i < length(text) - length(pattern) + 1; ++i)
        score[i] = computeLocalScore(infix(text, i, i + length(pattern)), pattern);

    return score;
}

int main()
{
    String<char> text = "This is an awesome tutorial to get to know SeqAn!";
    String<char> pattern = "tutorial";
    String<int> score = computeScore(text, pattern);

    for (unsigned i = 0; i < length(score); ++i)
        std::cout << score[i] << " ";
    std::cout << std::endl;

    return 0;
}

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.

The Role of References in SeqAn¶

Let us now have a closer look at the signature of computeScore().

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.

Therefore we change the signature of computeScore to:

String<int> computeScore(String<char> const & text, String<char> const & pattern)

Reading from right to left the function expects two references to const objects of type String of char.

Generic and Reusable Code¶

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:

template <typename TText, typename TPattern>
String<int> computeScore(TText const & text, TPattern const & pattern)

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.

From Object-Oriented Programming to SeqAn¶

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.

template <typename TText>
int computeLocalScore(TText const & subText, seqan::String<seqan::AminoAcid> const & pattern)
{
    int localScore = 0;
    for (unsigned i = 0; i < seqan::length(pattern); ++i)
        localScore += seqan::score(seqan::Blosum62(), subText[i], pattern[i]);

    return localScore;
}

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>, ...

struct SpecA;
struct SpecB;
struct SpecC;

template <typename TAlphabet, typename TSpec>
class String{};

template <typename TAlphabet, typename TSpec>
void myFunction(String<TAlphabet, TSpec> const &){}  // Variant (A)

template <typename TAlphabet>
void myFunction(String<TAlphabet, SpecB> const &){}  // Variant (B)

// ...

int main()
{
    String<char, SpecA> a;
    String<char, SpecB> b;
    String<char, SpecC> c;

    myFunction(a);            // calls (A)
    myFunction(b);            // calls (B)
    myFunction(c);            // calls (A)
}

Tags in SeqAn¶

Sometimes you will see something like this:

    globalAlignment(align, seqan::MyersHirschberg());

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.

template <typename T>
struct Tag{};

struct TagA_;
typedef Tag<TagA_> TagA;

struct TagB_;
typedef Tag<TagB_> TagB;

void myFunction(TagA const &){}  // (1)
void myFunction(TagB const &){}  // (2)

int main()
{
    myFunction(TagA());  // (3)
    myFunction(TagB());  // (4)
    return 0;
}

The code for the two different print() functions mentioned above could look like this:

#include <iostream>
#include <seqan/sequence.h>
#include <seqan/score.h>

template <typename TText, typename TSpec>
void print(TText const & text, TSpec const & /*tag*/)
{
    for (unsigned i = 0; i < seqan::length(text); ++i)
        std::cout << text[i] << " ";
    std::cout << std::endl;
}

struct MaxOnly {};

template <typename TText>
void print(TText const & score, MaxOnly const & /*tag*/)
{
    int maxScore = score[0];
    seqan::String<int> output;
    appendValue(output, 0);
    for (unsigned i = 1; i < seqan::length(score); ++i)
    {
        if (score[i] > maxScore)
        {
            maxScore = score[i];
            clear(output);
            resize(output, 1, i);
        }
        else if (score[i] == maxScore)
        {
            appendValue(output, i);
        }
    }

    for (unsigned i = 0; i < seqan::length(output); ++i)
        std::cout << output[i] << " ";
    std::cout << std::endl;
}

int main()
{
    return 0;
}

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.

The Final Result¶

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>

using namespace seqan;

template <typename TText>
int computeLocalScore(TText const & subText, String<AminoAcid> const & pattern)
{
    int localScore = 0;
    for (unsigned i = 0; i < length(pattern); ++i)
        localScore += score(Blosum62(), subText[i], pattern[i]);

    return localScore;
}

template <typename TText, typename TPattern>
int computeLocalScore(TText const & subText, TPattern const & pattern)
{
    int localScore = 0;
    for (unsigned i = 0; i < length(pattern); ++i)
        if (subText[i] == pattern[i])
            ++localScore;

    return localScore;
}

template <typename TText, typename TPattern>
String<int> computeScore(TText const & text, TPattern const & pattern)
{
    String<int> score;
    resize(score, length(text) - length(pattern) + 1, 0);

    for (unsigned i = 0; i < length(text) - length(pattern) + 1; ++i)
        score[i] = computeLocalScore(infix(text, i, i + length(pattern)), pattern);

    return score;
}

template <typename TText>
void print(TText const & text)
{
    std::cout << text << std::endl;
}

void print(String<int> const & text)
{
    for (unsigned i = 0; i < length(text); ++i)
        std::cout << text[i] << " ";
    std::cout << std::endl;
}

template <typename TText, typename TSpec>
void print(TText const & text, TSpec const & /*tag*/)
{
    print(text);
}

struct MaxOnly {};

template <typename TText>
void print(TText const & score, MaxOnly const & /*tag*/)
{
    int maxScore = score[0];
    String<int> output;
    appendValue(output, 0);
    for (unsigned i = 1; i < length(score); ++i)
    {
        if (score[i] > maxScore)
        {
            maxScore = score[i];
            clear(output);
            resize(output, 1, i);
        }
        else if (score[i] == maxScore)
            appendValue(output, i);
    }

    print(output);
}

struct GreaterZero {};

template <typename TText>
void print(TText const & score, GreaterZero const & /*tag*/)
{
    String<Pair<int> > output;
    for (unsigned i = 1; i < length(score); ++i)
        if (score[i] > 0)
            appendValue(output, Pair<int>(i, score[i]));

    for (unsigned i = 0; i < length(output); ++i)
        std::cout << "(" << output[i].i1 << "; " << output[i].i2 << ") ";
    std::cout << std::endl;
}

int main()
{
    String<char> text = "This is an awesome tutorial to get to now SeqAn!";
    String<char> pattern = "tutorial";
    String<int> score = computeScore(text, pattern);

    print(text);
    // > This is an awesome tutorial to get to now SeqAn!
    print(pattern);
    // > tutorial
    print(score);
    // > 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 8 0 1 0 0 0 0 2 0 1 0 0 1 0 3 0 1 1 0 0 0 0
    print(score, MaxOnly());
    // > 19
    print(score, GreaterZero());
    // > (2; 1) (5; 1) (12; 1) (17; 1) (19; 8) (21; 1) (26; 2) (28; 1) (31; 1) (33; 3) (35; 1) (36; 1)

    // And now for a protein pattern
    String<AminoAcid> protein = "tutorial";
    String<int> proteinScore = computeScore(text, protein);

    print(text);
    // > This is an awesome tutorial to get to now SeqAn!
    print(protein);
    // > TXTXRIAL
    print(proteinScore);
    // > 6 -9 -3 -6 -6 0 -9 -8 -7 -3 -9 -5 -8 -4 -5 -6 -6 1 -6 25 -7 2 -6 -6 -9 -6 -5 -7 1 -7 -5 -4 -6 2 -6 -3 -8 -9 -10 -4 -6 0 0 0 0 0 0 0
    print(proteinScore, MaxOnly());
    // > 19
    print(proteinScore, GreaterZero());
    // > (17; 1) (19; 25) (21; 2) (28; 1) (33; 2)

    return 0;
}

A First Working Example¶

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>

int main(int argc, char const ** argv)
{
    // Setup ArgumentParser.
    seqan::ArgumentParser parser("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::ParseResult res = 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)
        return res == seqan::ArgumentParser::PARSE_ERROR;

    // Extract option values and print them.
    unsigned period = 0;
    getOptionValue(period, parser, "period");
    bool toUppercase = isSet(parser, "uppercase");
    seqan::CharString text;
    getArgumentValue(text, parser, 0);

    std::cout << "period   \t" << period << '\n'
              << "uppercase\t" << toUppercase << '\n'
              << "text     \t" << text << '\n';

    return 0;
}

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

For example, we can already let the program generate an online help:

# modify_string -h
modify_string
=============

SYNOPSIS

DESCRIPTION
    -h, --help
          Displays this help message.
    -i, --period INT
          Period to use for the index.
    -U, --uppercase
          Select to-uppercase as operation.

VERSION
    modify_string version:
    Last update

While already informative, the help screen looks like there is something missing. For example, there is no synopsis, no version and no date of the last update given. We will fill this in later.

When we pass some parameters, the settings are printed:

# modify_string "This is a test." -i 1 -U
period     1
uppercase  1
text       This is a test.

When we try to use the --lowercase/-L option, we get an error. This is not surprising since we did not tell the argument parser about this option yet.

# modify_string "This is a test." -i 1 -L
modify_string: illegal option -- L

A Detailed Look¶

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::ArgumentParser parser("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.

    addArgument(parser, seqan::ArgParseArgument(
        seqan::ArgParseArgument::STRING, "TEXT"));

Then, we add options to the parser using addOption. We pass the parser and an ArgParseOption object.

    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.

Next, we parse the command line using parse.

    // Parse command line.
    seqan::ArgumentParser::ParseResult res = seqan::parse(parser, argc, argv);

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).

    if (res != seqan::ArgumentParser::PARSE_OK)
        return res == seqan::ArgumentParser::PARSE_ERROR;

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.

    // Extract option values and print them.
    unsigned period = 0;
    getOptionValue(period, parser, "period");
    bool toUppercase = isSet(parser, "uppercase");
    seqan::CharString text;
    getArgumentValue(text, parser, 0);

    std::cout << "period   \t" << period << '\n'
              << "uppercase\t" << toUppercase << '\n'
              << "text     \t" << text << '\n';

# program -a 1 -a 2 -a 3

If the option a is not a list then the occurence -a 3 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.

Using Default Values¶

Would it not be nice if we could specify a default value for --period, so it is 1 if not specified and simply each character is modified? We can do this by using the function setDefaultValue:

setDefaultValue(parser, "period", "1");

Note that we are giving the default value as a string. The ArgumentParser object will simply interpret it as if it was given on the command line. There, of course, each argument is a string.

Best Practice: Using Option Structs¶

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++.

The ModifyStringOptions struct’s looks as follows:

struct ModifyStringOptions
{
    unsigned period;
    bool toUppercase;
    bool toLowercase;
    seqan::CharString text;

    ModifyStringOptions() :
    period(1), toUppercase(false), toLowercase(false)
    {}
};

Click more... to see the whole updated program.

#include <iostream>

#include <seqan/arg_parse.h>

struct ModifyStringOptions
{
    unsigned period;
    bool toUppercase;
    bool toLowercase;
    seqan::CharString text;

    ModifyStringOptions() :
    period(1), toUppercase(false), toLowercase(false)
    {}
};

int main(int argc, char const ** argv)
{
    // Setup ArgumentParser.
    seqan::ArgumentParser parser("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::ParseResult res = 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)
        return res == seqan::ArgumentParser::PARSE_ERROR;

    // Extract option values and print them.
    ModifyStringOptions options;
    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';

    return 0;
}

Best Practice: Wrapping Parsing In Its Own Function¶

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>

struct ModifyStringOptions
{
    unsigned period;
    bool toUppercase;
    bool toLowercase;
    seqan::CharString text;

    ModifyStringOptions() :
    period(1), toUppercase(false), toLowercase(false)
    {}
};

seqan::ArgumentParser::ParseResult
parseCommandLine(ModifyStringOptions & options, int argc, char const ** argv)
{
    // Setup ArgumentParser.
    seqan::ArgumentParser parser("modify_string");

    // We require one argument.
    addArgument(parser, seqan::ArgParseArgument(
        seqan::ArgParseArgument::STRING, "TEXT"));

    // Define 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."));

    // Parse command line.
    seqan::ArgumentParser::ParseResult res = seqan::parse(parser, argc, argv);

    // Only extract  options if the program will continue after parseCommandLine()
    if (res != seqan::ArgumentParser::PARSE_OK)
        return res;

    // 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";
        return seqan::ArgumentParser::PARSE_ERROR;
    }

    return seqan::ArgumentParser::PARSE_OK;
}


int main(int argc, char const ** argv)
{
    // Parse the command line.
    ModifyStringOptions options;
    seqan::ArgumentParser::ParseResult res = 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)
        return res == 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';

    return 0;
}

Feature-Complete Example Program¶

The command line parsing part of our program is done now. Let us now add a function modifyText() 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>

struct ModifyStringOptions
{
    unsigned period;
    bool toUppercase;
    bool toLowercase;
    seqan::CharString text;

    ModifyStringOptions() :
    period(1), toUppercase(false), toLowercase(false)
    {}
};

seqan::ArgumentParser::ParseResult
parseCommandLine(ModifyStringOptions & options, int argc, char const ** argv)
{
    // Setup ArgumentParser.
    seqan::ArgumentParser parser("modify_string");

    // We require one argument.
    addArgument(parser, seqan::ArgParseArgument(
        seqan::ArgParseArgument::STRING, "TEXT"));

    // Define 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."));

    // Parse command line.
    seqan::ArgumentParser::ParseResult res = seqan::parse(parser, argc, argv);

    // Only extract  options if the program will continue after parseCommandLine()
    if (res != seqan::ArgumentParser::PARSE_OK)
        return res;

    // 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";
        return seqan::ArgumentParser::PARSE_ERROR;
    }

    return seqan::ArgumentParser::PARSE_OK;
}

seqan::CharString modifyString(seqan::CharString const & text,
                               ModifyStringOptions const & options)
{
    seqan::CharString result;

    if (options.toLowercase)
    {
        for (unsigned i = 0; i < length(text); ++i)
        {
            if (i % options.period == 0u)
                appendValue(result, tolower(text[i]));
            else
                appendValue(result, text[i]);
        }
    }
    else
    {
        for (unsigned i = 0; i < length(text); ++i)
        {
            if (i % options.period == 0u)
                appendValue(result, toupper(text[i]));
            else
                appendValue(result, text[i]);
        }
    }

    return result;
}

int main(int argc, char const ** argv)
{
    // Parse the command line.
    ModifyStringOptions options;
    seqan::ArgumentParser::ParseResult res = 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)
        return res == seqan::ArgumentParser::PARSE_ERROR;

    std::cout << modifyString(options.text, options) << '\n';

    return 0;
}

Setting Restrictions¶

One nice feature of the ArgumentParser is that it is able to perform some simple checks on the parameters. We can:

• check numbers for whether they are greater/smaller than some limits,
• mark options as being required, and
• setting lists of valid values for each option.

In this section, we will give some examples.

seqan::ArgumentParser parser("modify_string");
addOption(parser, seqan::ArgParseOption(
    "i", "integer-value", "An integer option",
    seqan::ArgParseArgument::INTEGER, "INT"));

setMinValue(parser, "i", "10");
setMaxValue(parser, "integer-value", "20");