 |
SeqAn3 3.4.0-rc.1
The Modern C++ library for sequence analysis.
|
|
SeqAn3 3.4.0-rc.1
The Modern C++ library for sequence analysis.
|
Learning Objective:
In this tutorial, you will learn how to construct an index and conduct searches.
| Difficulty | Moderate |
|---|---|
| Duration | 60 Minutes |
| Prerequisite tutorials | Ranges (Recommended) Pairwise Alignment (only for last assignment) |
| Recommended reading | FM-Index paper FM-Index on Wikipedia Bi-FM-Index paper |
Exact and approximate string matching is a common problem in bioinformatics, e.g. in read mapping. Usually, we want to search for a hit of one or many small sequences (queries) in a database consisting of one or more large sequences (references). Trivial approaches such as on-line search fail at this task for large data due to prohibitive run times.
The general solution is to use a data structure called index. An index is built over the reference and only needs to be computed once for a given data set. It is used to speed up the search in the reference and can be re-used by storing it to disk and loading it again without recomputation.
There are two groups of indices: hash tables and suffix-based indices. SeqAn implements the FM-Index and Bi-FM-Index as suffix-based indices. While the Bi-FM-Index needs more space (a factor of about 2), it allows faster searches.
Given an index, SeqAn will choose the best search strategy for you. Since very different algorithms may be selected internally depending on the configuration, it is advisable to do benchmarks with your application. A rule of thumb is to use the Bi-FM-Index when allowing more than 2 errors.
This tutorial will show you how to use the seqan3::fm_index and seqan3::bi_fm_index to create indices and how to search them efficiently using seqan3::search.
With this module you can:
The results of the search can be passed on to other modules, e.g. to create an alignment.
A reference is the data you want to search in, e.g. a genome or protein database.
A query is the data you want to search for in the reference, e.g. a read or an amino acid sequence.
An index is a data structure built over the reference that allows fast searches.
The full-text minute index (FM-Index) is an index that is similar to a suffix tree, but much smaller. It is used by most state-of-the-art read mappers and aligners. You can find more information on FM-Indicies in the original publication and on Wikipedia.
The bidirectional FM-Index (Bi-FM-Index) is an extension of the FM-Index that enables faster searches, especially when allowing multiple errors. But it uses almost twice the amount of memory the FM-Index uses. You can find more information on Bi-FM-Indicies here.
Constructing a (Bi-)FM-Index is very simple:
You can also index text collections (e.g. genomes with multiple chromosomes or protein databases):
The indices can also be stored and loaded from disk by using cereal.
Note that in contrast to the construction via a given text, the template cannot be deduced by the compiler when using the default constructor so you have to provide template arguments.
Create a seqan3::fm_index over the reference, store the index and load the index into a new seqan3::fm_index object. Print whether the indices are identical or differ.
Expected output:
Using an index, we can now conduct searches for a given query. In this part, we will learn how to search exactly, allow substitutions and indels, and how to configure what kind of results we want, e.g. all results vs. only the best result.
Exact search: Finds all locations of the query in the reference without any errors.
Approximate search: Finds all locations of the query in the reference with substitutions and indels within the confined set of maximal allowed errors.
Hit: A single result that identifies a specific location in the reference where a particular query can be found by either an exact or an approximate search.
We can search for all exact hits using seqan3::search:
You can also pass multiple queries at the same time:
The returned result is a lazy range over individual results, where each entry represents a specific location within the reference sequence for a particular query.
GCT in the text from assignment 1.
Expected output:
Up until now, we have seen that we can call the search with the query sequences and the index. In addition, we can provide a third parameter to provide a user defined search configuration. If we do not provide a user defined search configuration, the seqan3::search_cfg::default_configuration will be used, which triggers an exact search, finding all hits for a particular query. In the following, we will see how we can change the behaviour of the search algorithm by providing a user defined search configuration.
You can specify the error configuration for the approximate search using the seqan3::search_cfg::max_error_* configuration. Here you can use a combination of the following configuration elements to specify exactly how the errors can be distributed during the search:
Each of the configuration elements can be constructed with either an absolute number of errors or an error rate depending on the context. These are represented by the following types:
By combining the different error types using the |-operator, we give you full control over the error distribution. Thus, it is possible to set an upper limit of allowed errors but also to refine the error distribution by specifying the allowed errors for substitutions, insertions, or deletions.
For example, to search for either 1 insertion or 1 deletion you can write:
Here, we restrict the approximate search to only allow one error. This can be then either an insertion or a deletion but not both, since that would exceed the total error limit. This basically means that the error counts/rates do not have to sum up to the total of errors allowed:
In the above example, we allow 2 errors, which can be any combination of 2 substitutions, 1 insertion and 1 deletion. Defining only the total will set all error types to this value, i.e. if the total error is set to an error count of 2, any combination of 2 substitutions, 2 insertions and 2 deletions is allowed. On the other hand, when defining any of the error types but no total, the total will be set to the sum of all error types. For example, if we would not specify a total error of 1 in the first example above, the total error would be set to 2 automatically. Hence, the search will also include approximate hits containing one insertion and one deletion.
GCT in the text from assignment 1.
Allow up to 1 substitution and print all occurrences.
Expected output:
Besides the max error configuration, you can specify the scope of the search algorithm. This means that you can control which hits should be reported by the search.
To do so, you can use one of the following seqan3::search_cfg::hit_* configurations:
In contrast to the max error configuration, which allows a combination of the different error configuration objects, the hit configuration can only exist once within one search configuration. Trying to specify more than one hit configuration in one search configuration will fail at compile time with a static assertion. Sometimes the program you write requires to choose between different hit configurations depending on a user given program argument at runtime. To handle such cases you can also use the dynamic configuration seqan3::search_cfg::hit. This configuration object represents one of the four hit configurations mentioned previously and can be modified at runtime. The following snippet gives an example for this scenario:
Note that the same rules apply to both the dynamic and static hit configuration. That is, it can be added via the |-operator to the search configuration but cannot be combined with any other hit configuration.
A closer look at the strata configuration reveals that it is initialised with an additional parameter called the stratum. The stratum can be modified even after it was added to the search configuration like the following example demonstrates:
Here we introduced a new concept when working with the seqan3::configuration object, which is much like the access interface of a std::tuple. Concretely, it is possible to access the stored configuration using the get<cfg_type>(cfg) interface, where cfg_type is the name of the configuration type we would like to access. The get interface returns a reference to the stored object that is identified by the given name. If you try to access an object which does not exist within the search configuration, a static assert will be emitted at compile time such that no invalid code can be generated.
using seqan3::get; before we can call the get interface in order to allow the compiler to find the correct implementation of it based on the passed argument. This is related to how C++ resolves unqualified lookup of free functions in combination with function templates using an explicit template argument such as the get interface does.So, the open question remains what the stratum actually does. In the above example, if the best hit found by the search for a particular query had an edit distance of 1, the strata strategy would report all hits with up to an edit distance of 2. Since in this example the total error number is set to 2, all hits with 1 or 2 errors would be reported.
GCT in the text from assignment 1.seqan3::search_cfg::strata{1}).
Expected output:
When calling the search algorithm, a lazy range over seqan3::search_result objects is returned. Each result object represents a single hit. This means that merely calling the seqan3::search algorithm will do nothing except configure the search algorithm based on the given search configuration, query and index. Only when iterating over the lazy search result range, the actual search for every query is triggered. We have done this automatically in the previous examples when printing the result to the seqan3::debug_stream which then invokes a range based iteration over the returned range or by using the std::ranges::distance algorithm. However, in many cases, we want to access the specific positions and information stored in the seqan3::search_result object to proceed with our application. Since some information might be more compute-intensive than others, there is a way to control what the final search result object will contain.
The behaviour of the search algorithm is further controlled through the seqan3::search_cfg::output_* configurations. The following output configurations exists:
Similarly to the max error configurations, you can arbitrarily combine the configurations to customise the final output. For example, if you are only interested in the position of the hit within the reference sequence, you can use the seqan3::search_cfg::output_reference_begin_position configuration. Instead, if you need access to the index where the hit was found, you can use the seqan3::search_cfg::output_index_cursor configuration.
In the final example, we will extend our previous search examples to also compute the alignment of the found hits and their respective reference infixes. To do so, we recommend working through the Pairwise Alignment tutorial first.
GCT in the text from assignment 1.seqan3::search_cfg::hit_all_best.
Expected output:
1.9.8