Indices

A substring index is a datatype which allows to seek efficiently for all occurrences of a pattern in a string or a set of strings. Substring indices are very efficient for the exact string matching problem, i.e. finding all exact occurrences of a pattern in a text or a text collection. Instead of searching through the text in O(n) like online-search algorithms do, a substring index looks up the pattern in sublinear time o(n). Substring indices are full-text indices, i.e. they handle all substrings of a text in contrast to inverted files or signature files, which need word delimiters. SeqAn contains data structures to create, hold and use substring indices. Based on a unified concept, SeqAn offers the following concrete implementations defined as specializations of Index:

Specialization	Description
IndexEsa	Enhanced suffix array (Abouelhoda et al., 2004)
IndexWotd	Lazy suffix tree, constructed node-by-node on-demand (Giegerich et al., 2003)
IndexDfi	Lazy suffix tree for frequency based string mining (Weese, Schulz, 2008)
IndexQGram	q-gram index
PizzaChili	An adapter for the  Pizza & Chili index API

Suffix Tree Interface

The unified concept allows the first three indices ( IndexEsa, IndexWotd, IndexDfi) to be accessed just like a suffix tree independently of its concrete implementation. To access this (virtual) suffix tree SeqAn offers various iterators.

Depth-First Search

In SeqAn a suffix tree (see definition here) can be accessed with special suffix tree iterators, which differ in the way the tree nodes are traversed. For many sequence algorithms it is neccessary to do 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 (Fig.1) halts in a node when visiting it for the first time whereas a postorder traversal (Fig.2) halts when visiting a node for the last time. The following two figures give an example in which order the tree nodes are visited.

Figure 1: Preorder DFS	Figure 2: Postorder DFS

There are currently 2 iterators in SeqAn supporting a DFS search:

Iterator	Preorder	Postorder
BottomUp Iterator	-	+
TopDownHistory Iterator	+	+

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.

Example

We want to construct the suffix tree of the string "abracadabra" and output the substrings represented by tree nodes in preorder dfs. All index related data structures are defined in seqan/index.h, so we have to include this header (seqan/sequence.h is included implicitly).

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

using namespace seqan;

Next we create the string "abracadabra" and an index specialized with the type of this string. The string can be given to the index constructor.

int main ()
{
    String<char> myString = "abracadabra";
    typedef Index< String<char> > TMyIndex;
    TMyIndex myIndex(myString);

The Iterator metafunction expects two arguments, the type of the container to be iterated and a specialization tag, see the VSTree Iterator hierarchy. In this example we chose a TopDownHistory Iterator whose signature in the second template argument is TopDown< ParentLinks<Preorder> >. Each suffix tree iterator constructor expects at least the index object and optional problem-dependent parameters. As all DFS suffix tree iterators implement the Iterator concept, they can be used via goNext, atEnd, etc. The string that represents the node the iterator points to is returned by representative.

    Iterator< TMyIndex, TopDown< ParentLinks<Preorder> > >::Type myIterator(myIndex);
    while (!atEnd(myIterator))
    {
        std::cout << representative(myIterator) << std::endl;
        ++myIterator;
    }
    return 0;
}

Program output:

a
abra
abracadabra
acadabra
adabra
bra
bracadabra
cadabra
dabra
ra
racadabra

Note: A relaxed suffix tree (see definition) 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.

Assignments

Task 1

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.

Difficulty

2

Solution

can be found here

Task 2

Write a program that outputs all maximal unique matches (MUMs) between "CDFGHC" and "CDEFGAHC".

Difficulty

2

Solution

can be found here

Top-Down Iteration

For index based pattern search or algorithms traversing only the upper parts of the suffix tree the TopDown Iterator or TopDownHistory Iterator is the best solution. Both provide the functions goDown and goRight to go down to the first child node or go to the next sibling. The TopDownHistory Iterator additionally provides goUp to go back to the parent node. The child nodes in IndexEsa indices are lexicographically sorted from first to last. For IndexWotd and IndexDfi indices this holds for all children except the first.

Example

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. First we create an index of the text "How many wood would a woodchuck chuck."

int main ()
{
    typedef Index<CharString> TIndex;
    TIndex index("How many wood would a woodchuck chuck.");
    Iterator< TIndex, TopDown<> >::Type it(index);

The main search can then 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.

    CharString pattern = "wood";
    while (repLength(it) < length(pattern))
    {
        // go down edge starting with the next pattern character
        if (!goDown(it, pattern[repLength(it)])) return 0;
        unsigned endPos = _min(repLength(it), length(pattern));
        // compare remaining edge characters with pattern
        std::cout << representative(it) << std::endl;
        if (infix(representative(it), parentRepLength(it) + 1, endPos) !=
            infix(pattern, parentRepLength(it) + 1, endPos)) return 0;
    }

If all pattern characters could successfully be compared we end in the topmost node pattern is a prefix of. Thus, the suffixes represented by this node are the occurrences of our pattern.

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

Program output:

w
wo
wood
9
22

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

    if (goDown(it, "wood"))
        for (unsigned i = 0; i < length(getOccurrences(it)); ++i)
            std::cout << getOccurrences(it)[i] << std::endl;
    return 0;
}

9
22

Assignments

Task 3

Write a program that iterates over all nodes of the suffix tree of the string "tobeornottobe" in preorder DFS. Use goDown, goRight and goUp to iterate instead of goNext or the operator++. Output the representatives.

Difficulty

4

Solution

can be found here

Task 4

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.

Difficulty

5

Solution

can be found here

Access Suffix Tree Nodes

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, so we concentrate on the most important ones.

Function	Description
representative	returns the substring that represents the current node, i.e. the concatenation of substrings on the path from the root to the current node
getOccurrence	returns a position where the representative occurs in the text
getOccurrences	returns a string of all positions where the representative occurs in the text
isRightTerminal	tests if the representative is a suffix in the text (corresponds to the shaded nodes in the suffix tree figures)
isLeaf	tests if the current node is a tree leaf
parentEdgeLabel	returns the substring that represents the edge from the current node to its parent (only TopDownHistory Iterator)

Note: There is a difference between the functions isLeaf and isRightTerminal. In a relaxed suffix tree (see definition) 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.

Property Maps

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

int main ()
{
    String<char> myString = "abracadabra";
    typedef Index< String<char>, IndexWotd<> > TMyIndex;
    TMyIndex myIndex(myString);
    String<int> propMap;

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.

    Iterator< TMyIndex, TopDown< ParentLinks<Preorder> > >::Type myIterator(myIndex);
    int depth;
    while (!atEnd(myIterator))
    {
        if (isRoot(myIterator))
            depth = 0;
        else
            depth = getProperty(propMap, nodeUp(myIterator)) + 1;
        resizeVertexMap(myIndex, propMap);
        assignProperty(propMap, value(myIterator), depth);
        ++myIterator;
    }

At the end we again iterate over all nodes and output the calculated node depth.

    goBegin(myIterator);
    while (!atEnd(myIterator))
    {
        std::cout << getProperty(propMap, value(myIterator)) << '\t' << representative(myIterator) << std::endl;
        ++myIterator;
    }
    return 0;
}

Program output:

0
1       a
2       abra
3       abracadabra
2       acadabra
2       adabra
1       bra
2       bracadabra
1       cadabra
1       dabra
1       ra
2       racadabra

Hint

In SeqAn there is already a function nodeDepth defined to return the node depth.

Additional iterators

By now, we know the following iterators (n=text size, σ=alphabet size, d=tree depth):

Iterator specialization	Description	Space	Index tables
BottomUp Iterator	postorder dfs	O(d)	SA, LCP
TopDown Iterator	can go down and go right	O(1)	SA, Lcp, Childtab
TopDownHistory Iterator	can also go up, preorder/postorder dfs	O(d)	SA, Lcp, Childtab

Besides the iterators described above, there are some application-specific iterators in SeqAn:

Iterator specialization	Description	Space	Index tables
MaxRepeats Iterator	maximal repeats	O(n)	SA, Lcp, Bwt
SuperMaxRepeats Iterator	supermaximal repeats	O(d+σ)	SA, Lcp, Childtab, Bwt
SuperMaxRepeatsFast Iterator	supermaximal repeats (optimized for enh. suffix arrays)	O(σ)	SA, Lcp, Bwt
MUMs Iterator	maximal unique matches	O(d)	SA, Lcp, Bwt
MultiMEMs Iterator	multiple maximal exact matches (w.i.p.)	O(n)	SA, Lcp, Bwt

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 s1, ..., sm a MultiMEM (multiple maximal exact match) is a substring r that occurs in each sequence si 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:

Example

Maximal Unique Matches

Supermaximal Repeats

Maximal Repeats

q-gram Index

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 HowTo), 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:

Specialization	Modifiable*	Number of Gaps
UngappedShape	-	0
SimpleShape	+	0
OneGappedShape	+	0/1
GappedShape	-	any
GenericShape	+	any

• - fixed at compile time, + can be changed at runtime

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=((3*4+3)*4+1. 4 is the alphabet size in this example (see ValueSize).

The q-gram index offers different function 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 as the ESA index and in the same fibre (Fibre_SA). If only the number of q-grams per sequence are needed the QGram_Counts and QGram_CountsDir 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 QGram_SA, QGram_Counts fibres), called buckets. The begin position of bucket i is stored in directory fibres (QGram_Dir, QGram_CountsDir) 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 OpenAddressing 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 (QGram_BucketMap).

For more details on q-gram index fibres see the HowTo or QGram Index Fibres.

Example

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

int main ()
{
    typedef Index<DnaString, IndexQGram< UngappedShape<3> > > TIndex;
    TIndex index("CATGATTACATA");

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.

    hash(indexShape(index), "CAT");
    for (unsigned i = 0; i < length(getOccurrences(index, indexShape(index))); ++i)
        std::cout << getOccurrences(index, indexShape(index))[i] << std::endl;
    return 0;
}

Program output:

0
8

Assignments

Task 5

Write a program that outputs all occurrences of the gapped q-gram "AT-A" in "CATGATTACATA".

Difficulty

3

Solution

can be found here

Task 6

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 OpenAddressing index and output the directory sizes (QGram_Dir, QGram_CountsDir fibre).

Difficulty

5

Hint

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

Solution

can be found here

Handling Multiple Sequences

The previous sections briefly described how an index of a set of strings can be instantiated. Instead of creating an Index of a String you create one of a StringSet. A character position of this string set can be one of the following:

1. A local position (default), i.e. Pair (seqNo, seqOfs) where seqNo identifies the string within the stringset and the seqOfs identifies the position within this string.
2. A global position, i.e. single integer value between 0 and the sum of string lengths minus 1 (global position). This integer is the position in the gapless concatenation of all strings in the StringSet to a single string.

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 getOccurrence or 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.

Submit a comment

If you found a mistake, or have suggestions about an improvement of this page press: submit your comment