When searching databases of nucleotide or protein sequences, finding a local alignment of two sequences is one of the main tasks. Since the sizes of available databases grow constantly, the efficiency of retrieval methods becomes the critical issue. The sequence retrieval relies on finding sequences in the database which align best with the query sequence. However, an optimal alignment can be found in quadratic time (by use of dynamic programming) while this is infeasible when dealing with large databases. The existing solutions use fast heuristic methods (like BLAST, FASTA) which produce only an uncontrolled approximation of the best alignment and even do not provide any information about the alignment approximation error. In this paper we propose an approach of exact and approximate indexing using several metric access methods (MAMs) in combination with the TriGen algorithm, in order to reduce the number of alignments (distance computations) needed. The experimental results have shown that a straightforward adoption of MAMs to sequence retrieval cannot outperform the specialized heuristic algorithms (at least at the moment). On the other side, the results show MAMs could provide a basis for specialized access methods capable of precision/efficiency trade-off control.
When managing biological sequence information, it is important to realize that organisms linked by common descent use similar genes to secure their biological functions. Therefore, it makes sense to store known genetic sequences in a database and search for similarities among them. Because these databases grow quickly in time and due to expensive comparison of sequences, there is a strong need for efficient methods capable of handling such amounts of data.
One way to handle the increasing number of sequences is indexing, which
is unavoidable in today’s situation of exponentially growing databases to be
searched (Figure 1). There have appeared indexing applications in this area
already, however, they usually tried to split sequences into q-grams (substrings
of length q) [
DNA molecules consist of a linear (unbranched) string of nucleotides (conventionally labeled A, C, G, T) which can be transcribed to RNA and later translated to proteins. A protein molecule is a linear chain of amino acids. There are 20 amino acids, each encoded by a triplet of nucleotides (called codon). The assignment between codons and amino acids is determined by an evolutionarily fixed code table – the genetic code. Every protein has some biological activity which is derived from its three-dimensional structure. Similar amino acid sequences tend to have similar three-dimensional structure and therefore similar function1. When a protein sequence is determined (nowadays it is usually derived from an experimentally obtained sequence of the corresponding gene), its function is usually unknown. One should try to find similar proteins in the database of already known protein sequences (even in proteins of different species) to find out possible purpose of the newly sequenced protein or at least to get a clue of it. Therefore, it is appropriate to search through as large amount of data as possible to increase the probability of finding similar protein. 1.2
Using nucleotide databases is not as common as using protein databases2. They are searched only for similarities among sequences of one species (unlike protein databases) which follows from the aims why they are analyzed, i.e. – finding similarities in parts of nucleotide sequences which are not transcribed to mRNA (non-coding sequences) – checking whether someone else has already sequenced given segment of DNA – checking whether given segment was sequenced incorrectly 1 Moreover more extensive sequence similarity can be viewed as evidence of common ancestry, and therefore as a basis for reconstructing phylogenetic history of organisms and their genes. 2 which is the reason why we have tested our method for protein sequences only, but it can be used for nucleotide sequences without any change at all.
As comes out from previous, the role of databases in bioinformatics is of major
importance. Nowadays, there exist three most prominent databases of sequences
– American GenBank [
Besides these databases, there exist several other (mostly protein) databases.
In our experiments we used the Swiss-Prot [
GenBank growth
Swiss-Prot growth 5 6 5 5 5 4 5 3 5 2 5 1 5 0
Base pairs of DNA (bil ions) Sequences (mil ions) – analysis of the dot matrix – dynamic programming solution – heuristic methods using q-grams
All these methods compute a kind of biological similarity between two sequences. The similarity of two sequences is defined by an optimal alignment of them. An alignment of two3 sequences is every such correspondence between letters in both (sub)sequences (one written under the other), which can be achieved by inserting gaps into either sequence. Naturally, among all possible alignments we are interested in an optimal one – which is the alignment with best realvalue score (interpreted as sequence similarity). A particular method of scoring 3 There can be also multiple sequences mutually aligned for some purposes, but this is out of scope of this paper. an alignment changes with method used (see next subsections), however, it is always a sum over the scores between aligned letters/gaps.
Consequently, the similarity assessment (score of alignment) can be used in retrieval of sequences from a database similar to a query sequence. 2.1
A dot matrix (dotplot) analysis is a method for finding possible alignments of two sequences. In the dot matrix analysis we mark each of the axes with one amino acid sequence and then mark each position in the matrix where the two sequences match (have the same letters on given positions). Matching segments can then be seen as diagonal rows of dots going from left (see Figure 2).4
This algorithm can be improved by use of sliding windows of given length5 (see Figure 2a,b). The window is successively placed to every position, while the window on every position is checked for a diagonal occurrence. If the number of matching letters in both sequences is above a fixed threshold (possibly according to a scoring matrix, see below), then a dot is placed on the position of the upper left corner of the window. Subsequently, short diagonals in the resulting plot can be filtered out for better transparency(Figure 2c).
Afterwards, the dot matrix is visually checked for diagonals long enough to be used as candidates for even longer alignment. Pairs of sequences showing sufficient similarity are then further inspected with some other techniques (e.g. tools using dynamic programming).
IR ! !! PIQITE !!!!!! SLAZ !!! ! !! TTEPFP !! ! !!! ! A (a)
IPIRTIQE ! !!!!!!!!! !! STETPLAFPZ !!!!!!! !!!!!!! !!!!! A
(b)
The dot matrix method is appropriate for visualisation and manual analysis of partial alignments of two sequences, but it does not give us the optimal alignment – especially if an alignment with gaps is our goal. To solve this problem, there exist algorithms based on dynamic programming which can compute optimal alignment with gaps in O(mn) (where m, n are sequence lengths) time. 4 This method can also be used to finding reverse matchings (diagonals go from right to left in that case). 5 When comparing nucleotide sequences, longer windows are used. With minor changes one can use the same algorithm for finding both, global and local alignment. A global alignment of sequences is such an alignment where the whole sequences must be aligned, whereas a local alignment methods look for best alignment for any subsequences of the input sequences.
First of all we need to define a distance to measure similarity between two sequences (possibly by use of an alignment). The simplest is the Hamming distance (HD) which is defined on two strings of equal lengths, and computes the number of positions for which the corresponding symbols are different. Because HD is a metric, it is often used when indexing sequences split to q-grams.
The Hamming distance is not suitable because it is sensitive to shifts in sequences. Therefore, a better measure is the Levenshtein (or edit) distance given by the minimum number of operations needed to transform one sequence into the other, where an operation is an insertion, deletion, or replacement of a single letter. The Levenshtein distance can be further extended with operation weights. Every replacement, insertion or deletion can be penalized by a constant value. Such an extended edit distance is called weighted edit distance.
However, when considering DNA or protein sequences, even weighted edit
distance is not suitable. The reason is that probability of mutations from one amino
acid into another (appearing in an alignment) is dependent on the two
particular amino acids involved (the same holds for nucleotides but this is not so often
used, since there are only four of them). To deal with this fact, a scoring matrix
is used which is a 20 × 20 square matrix where on position [i, j] is the so-called
weight of mutation of amino acid i into amino acid j (if we use identity matrix
we get an ordinary edit distance). There are many different scoring matrices used
for different purposes (e.g. PAM matrices [
si,j−1 + σ si,j + δ(ai, bj ) (1) (2) where a and b represent the sequences to be aligned, σ is a score for gaps and δ is scoring matrix. Since [i, j] contains the score of global alignment of the i-long prefix of a and j-long prefix of b, the cell [|a|, |b|] contains the alignment score for the whole sequences a, b.
Local Alignment Measures. Since we are more interested in matching
subsequences, the local alignment is applied more often in computing similarity
between biological sequences6. Local alignment of sequences a and b is finding a
subsequence S(a) of sequence a and a subsequence S(b) of sequence b, such that
global alignment of S(a) and S(b) provides the highest score. The algorithm for
local alignment can be seen as a modification of Needleman-Wunsch and was
published in 1981 by Smith and Waterman [
Finally, there is usually another modification which enables differentiating between the opening and extending gap. Extending a gap is considerable less penalized than opening a gap (i.e. two single gaps in alignment are more penalized than a single two-letter gap). 2.2
The Smith-Watterman (SW) algorithm gives the optimal solution to local
alignment but it is computationally expensive. Hence, there have been developed
cheaper heuristics which approximate the optimal local alignment. The first
wide-spread method of this type was the FASTA [
Such alignments are called high scoring pairs (HSP). 7. All HSPs with scores below a given threshold are excluded. 8. The scores of non-filtered sequences are refined by the classic Smith-Watterman algorithm.
Statistical Relevance. Because of the heuristic nature of BLAST, it is important to have a method saying how good the found hits are (however, there is no way to determine the retrieval error, say precision/recall in IR terminology). We are also interested in determining whether the retrieved alignments are some coincidence, or whether they really refer to a biological relation. The local alignments are statistically well understood, so we are able to say how relevant the retrieved sequences are. This statistic evidence fits primarily to ungapped local alignment, but it has been shown (mostly by computational experiments) that it applies to gapped alignments, too.
We are interested whether the resulting score (either Smith-Watterman (SW)
or any other) has any statistical significance (otherwise the score alone is just a
number), hence, we have to take the score distribution into account, considering
the entire sequence space and a particular scoring matrix. We can derive an
expected number of sequences of lengths m and n with score at least S as
E = Kmne−λS
(3)
where K and λ are characteristics of the SW score distribution. The E formula
is called the E-value for the score S [
The E-value as defined above does not take into consideration dealing with more than two sequences, however, we should additionally take into account the fact that probability of finding a sequence with given score depends also on the number of sequences in the database or/and the query length. In order to prevent the effect that finding a short sequence with given E-value has the same probability as finding a longer sequence with the same E-value, we should multiply the E-value by the total number of residues (letters) in the database and obtain a new E-value (used by FASTA). Alternatively, instead of multiplying the E-value by database size, we could adjust the E-value by distinguishing different sequence lengths, since a query is more likely to be related to a long sequence than to a shorter one (used by BLAST). 7 This applies to BLAST2 - previous versions of BLAST did not connect sequences on diagonals (and therefore used higher value for cut-off score in step 4).
Since a particular method of sequence alignment can be viewed as a non-metric
distance/similarity measure, we had an idea to turn it into a distance metric
δ, which satisfies the metric properties (reflexivity, non-negativity, symmetry
and triangle inequality). Such a metric could be then utilized by various metric
access methods (MAMs) which have been designed to quickly search in databases
modeled in metric spaces [
First, we have to decide which sequence similarity measure should be used. The simplest solution would be the Smith-Waterman (SW) score. However, the problem with SW is that it is similarity measure, not a distance (greater scores mean higher similarities). Another problem is that SW lacks statistical significance. Therefore, the use of E-value is much better, since it is a distance and it is also widely used in ranking the results retrieved by FASTA or BLAST. Due to this choice, we can also compare our approach to the other methods.
If we want to use E-value for metric indexing, we need an equivalent metric. As defined, the original E-value satisfies just the non-negativity property, but it can be easily modified to satisfy also reflexivity and symmetry. To enforce reflexivity, the definition of E-value can be modified in a way that identical sequences have a zero E-value. This is not a problem because two identical sequences are also most similar in mathematical and biological meaning (no other sequence can be more similar). To satisfy the symmetry, we have to accomplish a more substantial change. As mentioned in Section 2.2, the BLAST’s modification of E-value takes into account the total number N of residues (letters) in database and the length n of inspected database sequence, as follows:
E = Kmne−λS N/n (4)
If we change the role of query and database sequence, we get different values. So we changed E-value to the following form:
E = Kmne−λS N/max(m, n) (5) where m is the length of query sequence. This is a minor change in statistical relevance if we realize that an average length of query is very similar to the average length of database sequences. As mentioned before, e.g. FASTA does not take sequence lengths into account at all.
TriGen Algorithm. Up to now, we have turned the E-value into a
semimetric (reflexive, non-negative, symmetric distance). The last step is to enforce
the triangle inequality, which is the hardest part. For this purpose we can use
the TriGen algorithm [
The LAESA method [
Whenever a query is to be processed, the query object is projected into the pivot space the same way as if it would be inserted. Then the distance matrix is sequentially searched and all the vectors which do not overlap the query region in the pivot space are filtered out. The remaining candidate objects (corresponding to the non-filtered vectors) are subsequently filtered in the original metric space.
The LAESA method is very powerful in its pruning effectiveness, however, due to expensive selection of pivots and due to the sequential processing of distance matrix its usage in dynamic database environments is limited. 3.3
A typical tree-based MAM designed for database environments is the M-tree [
A range query in M-tree is processed by processing just the nodes the region ball of which overlaps the query ball. A kNN query processing is similar, however, the radius of query ball is unknown at the beginning, so it must be heuristically updated during query evaluation.
As a combination of LAESA and M-tree, the PM-tree [
Slim-down Algorithm. A particular hierarchy of M-tree’s (or PM-tree’s)
nested data regions can be far from optimal, which means the volumes of
regions can be very large, and they can also overlap significantly. This leads to
poor search efficiency, since the larger volumes/overlaps the greater
probability that a data region and the query region will overlap. To prevent such poor
hierarchies, the generalized Slim-Down algorithm [
As dataset we used random subset of the Swiss-Prot database of size 3000 with total number of 1041027 amino acids. Another random 100 hundred sequences have been chosen as query sequences.All of the sequences were of maximal length 1000 which is doesn’t cause any problem when we realize that there are only 9191 longer sequences out of 252616 in whole Swiss-prot which makes 3% (average sequence length is 365 in whole Swiss-prot and 335 in the reduced variant). These longer sequences could then be treated in special way since they are just small part of the whole.
We don’t show time comparation in our tests because we do not have effecitive implementation of Smith-Waterman yet, being the crucial component of the running time. But our method could be easily compared to SSEARCH (part of FASTA package) when we realize that SSEARCH is equivalent to sequence scan. To be able to compare index based methods with BLAST we distinguish number of distance computations from computational costs. We defined computational costs here as number of comparing two letters. Therefore computational cost of tree based methods are averaged as number of distance computations multiplied by the average size of distance matrix for Smith-Watterman which is 321 ∗ 321 = 103041. For possibility of comparing index-based methods with sequentional scan, we show number of distance computations too. Computational costs of BLAST can be devided from description in section 2.2. Average number of neighbouring words for our dataset of query sequences is 54 and takes 81784 operations (matching letters) to get them 8 , which gives us 245352 computational operations. This means that there was built search tree with 54 items. Such a tree has aproximatly 6 levels. When searching for seeds for the high scoring pairs, there have to be done 6 comparations with every position of the database. Since every comparation of word with database costs 3 operations, the computation costs of BLAST in average are 6 ∗ 1041000 ∗ 3 = 1873800 operations. Number of operations done while finding neighbouring words and finding high scoring pairs is insignificant in comparation with this number so we do not take it into account. Therefore we state that in average case, computational costs of BLAST are 18738000 operations9.
Four indexing methods were tested - M-tree, PM-tree, slimmed PM-tree and LAESA, each of them using the same set of distance modificators 10 genereted by the TriGen algorithm. As can be seen in Figure 5a, when the zero error tolerance is used, weight of the modificator causes that number of distance computations is almost equivalent to sequence scan. When the weight is too big, it makes trinagle inequvality hold but for the price of increasing intrinsic dimension. However M-tree performes slightly better then sequence scan which means that when searching in the tree, not all of the nodes have been inspected (because of inner nodes, number of objects in the tree exceeds 3000). On the other side, PM-tree and slimmed PM-tree show worse result since there are aditional computations to pivots (mapping the query). The most similar to sequential scan is LAESA method, which uses constant number of distance computations independently on range of the query.
The situation slightly changes when we allow small error (Figure 5b). This causes that number of distance computations decrease about six percent compared to zero error tolerance. Here, PM-tree and slimmed PM tree outperform M-tree and the difference is about 1%. On the other hand LAESA showed just 8 found out empirically 9 We do not consider possible filtering of nonsignificant segments. 10 fractional power modificators were used
Error tolerance 0
Error tolerance - 0.002 2 slight improvement. But in both cases BLAST method is evidently more effective since effctivity of the index is almost sequence scan even if small error is allowed.
Why PM-tree and slimmed PM-tree behave better when allowing some error? Answer to this question can be seen on Figure 5c which shows on range query of E-value five the relation between declared TriGen error tolerance and real error experienced in test. Here we can see, that PM-tree real error growth more quickly than the error of M-tree and moreover, we can see that slope of those lines are almost inverse, which means that the distance computations gain of PM-tree is counterbalanced by the error. From these two graphs can also be seen that for real error 50%, there still have to be done approximately 1200 distance computations which is about 120000000 computational operations. That means that if BLAST would be such a bad hauristic that it would have just 50% successfulness, it still would be noticeable more effective. 5
In this paper, we have have tested suitability of metric access indexing methods for indexing protein sequences. It has been shown that these method are not applicable to sequence alignment problem without their modification. This is primarily because of quality of the data to be indexed and the distance function which is used to define similarity between them. This distance function is highly non-metric which demands strong modifications to it, to make it metric. This modification distorts distances in a way that strongly increases intrinsic dimension of the data and therefore the efficiency of search is almost the same as efficiency of sequential scan. But against sequentional scan it has that advantage that precission can be defined and thus traded off for efficiency. This learned facts can aim next resarch to several areas. To name a few: – Examining TriGen modificators and finding such modificators, which would minimize real error while distributing objects (i.e. sequences) in the space in a way which will be appropriate for indexing methods (i.e. descreasing intrinsic dimension). – Modifying the search structures. For example examining possiblities of cutting sequences to q-grams but being able to define arised error (caused by splitting and thus losing information included in the whole sequence) and (optimaly) minimalize it. – Modifying computing of Smith-Waterman local alignment. The idea is to change computing so that it will be faster and resulting scores won’t violate properites of metric so much, as they do now (for example by using borders to limit the computational space in the distance matrix).
This research has been partially supported by GACˇ R grant 201/05/P036 provided by the Czech Science Foundation. We would like to thank prof. Fatima Cvrˇckov´a (Department of Plant Physiology, Faculty of Science, Charles University in Prague) for helping to get in touch with biologist point of view to the problem and contributing to the paper with helpful comments.