String Indices

Learning Objective
You will get an overview of the different kinds of indices in SeqAn and how they are used.
Difficulty
Average
Duration
1 h
Prerequisites
Sequences

Indices in SeqAn

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

You can find the following indices in SeqAn:

IndexSa
Suffix Array [MM93]
IndexEsa
Extended Suffix Array [AKO04]
IndexWotd
Lazy suffix tree [GKS03]
IndexDfi
Deferred Frequency Index [WS08]
IndexQGram
Q-gram index (see here for more details)
FMIndex
Full-text minute index (see the FMIndex for more details) [FM01]

Index Construction

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

# export TMPDIR=/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.

    String<Dna5> genome = "ACGTACGTACGTN";
    Index<String<Dna5>, IndexEsa<> > esaIndex(genome);

In contrast, the next code snipped creates a FM index over a set of amino acid sequences:

    StringSet<String<AminoAcid> > protein;
    appendValue(protein, "VXLAGZ");
    appendValue(protein, "GKTVXL");
    appendValue(protein, "XLZ");

    Index<StringSet<String<AminoAcid> >, FMIndex<> > fmIndex(protein);

Assignment 1

Type
Review
Objective

Copy the code below and

  1. change it to build an IndexEsa over a string of type Dna,
  2. add an IndexEsa over a StringSet of Strings of type Dna.
#include <seqan/sequence.h>
#include <seqan/index.h>

using namespace seqan;

int main()
{
    String<char> text = "This is the first example";
    Index<String<char>, FMIndex<> > index(text);

    return 0;
}
Solution
#include <seqan/sequence.h>
#include <seqan/index.h>

using namespace seqan;

int main()
{
    // One possible solution to the first sub assignment
    String<Dna> text = "ACGTTTGACAGCT";
    Index<String<Dna>, IndexEsa<> > index(text);

    // One possible solution to the second sub assignment
    StringSet<String<Dna> > stringSet;
    appendValue(stringSet, "ACGTCATCAT");
    appendValue(stringSet, "ACTTTG");
    appendValue(stringSet, "CACCCCCCTATTT");

    Index<StringSet<String<Dna> >, IndexEsa<> > indexSet(stringSet);

    return 0;
}

Assignment 2

Type
Application
Objective
Write a small program that prints the locations of all occurrences of "TATAA" in "TTATTAAGCGTATAGCCCTATAAATATAA".
Hints
Use the find function as the conditional instruction of a <tt>while</tt> loop.
Solution
#include <seqan/sequence.h>
#include <seqan/index.h>

using namespace seqan;

int main()
{
    String<Dna5> genome = "TTATTAAGCGTATAGCCCTATAAATATAA";
    Index<String<Dna5>, IndexEsa<> > esaIndex(genome);
    Finder<Index<String<Dna5>, IndexEsa<> > > esaFinder(esaIndex);

    while (find(esaFinder, "TATAA"))
    {
        std::cout << position(esaFinder) << std::endl;
    }

    return 0;
}

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.

Handling Multiple Sequences (StringSets)

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:

  1. 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.
  2. 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.

Storing and Loading

Storing and loading an index can be done with:

    const char * tempFileName = SEQAN_TEMP_FILENAME();
    save(index, tempFileName);

or

    open(index, tempFileName);

If you have built your q-gram index with variable shapes (i.e. SimpleShape GenericShape), 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.

    indexRequire(index, FibreSA());

For the IndexEsa index you could do:

    indexRequire(esaIndex, EsaSA());
    indexRequire(esaIndex, EsaLcp());
    indexRequire(esaIndex, EsaChildtab());  // for TopDown iterators
    indexRequire(esaIndex, EsaBwt());       // for (Super-)MaxRepeats iterators

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:

    Index<String<Dna, External<> >, IndexEsa<> > external_index;
    open(external_index, tempFileName);

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.

Reducing the memory consumption

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.

SA Fibre

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<>
struct SAValue<String<Dna> >
{
    typedef unsigned Type;
};

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

template<>
struct SAValue<StringSet<String<Dna> > >
{
    typedef Pair<unsigned, unsigned> Type;
};

The first type of the pair is used as the type for the index of a string in the string set. So if you only have a few strings you could save even more memory like this.

typedef Pair<unsigned char, unsigned> Type;

FMIndex Fibres

The size of a generalized FMIndex depends also on the total number of characters in a StringSet (see lengthSum). This trait can be configured via the FMIndexConfig object. For more information, see the FMIndex section.

    typedef FMIndexConfig<void, unsigned> TConfig;
    Index<StringSet<String<Dna> >, FMIndex<void, TConfig> > configIndex(set);

Other Index Fibres

See Accessing Index Fibres for more information.