@ewha-bio:11 JSONTXT

A Heuristic Algorithm to Find All Normalized Local Alignments Above Threshold. Local alignment is an important task in molecular biology to see if two sequences contain regions that are similar. The most popular approach to local alignment is the use of dynamic programming due to Smith and Waterman, but the alignment reported by the Smith-Waterman algorithm has some undesirable properties. The recent approach to fix these problems is to use the notion of normalized scores for local alignments by Arslan, Egecioglu and Pevzner. In this paper we consider the problem of finding all local alignments whose normalized scores are above a given threshold, and present a fast heuristic algo­ rithm. Our algorithm is 180-330 times faster than Arslan et al.’s for sequences of length about 120 kbp and about 40-50 times faster for sequences of length about 30 kbp. Sequence alignment is a fundamental task in molecular biology to see if two sequences are related. The local alignment problem for two sequences is to find a pair of similar regions, one from each sequence. In many biological applications local alignment is more useful than global alignment because two sequences may not be similar as a whole, but they can contain small regions that are very similar. A most popular approach to local alignment is the use of dynamic programming due to Smith and Waterman (Smith and Waterman, 1981; Waterman et al., 1995). The Smith- Waterman algorithm finds a pair of regions whose alignment score is the highest (which is called the optimal local alignment). However, the alignment reported by the Smith-Waterman algorithm has some undesirable properties. The alignment may include regions whose similarity is very poor (Arslan et al., 2001; Zhang et al., 1999). Also, the Smith-Waterman algorithm does not consider the lengths of local alignments in computing scores, and therefore a biologically important but short alignment may not be detected(Arslan and Pevzner et al., 2001). There have been several approaches to fix the problems of the Smith-Waterman algorithm (Goad and Kanehisa, 1982; Sellers, 1984; Zhang eta/., 1998; Zhang etal., 1999; Arslan and Pevzner et a/., 2001). For example, Zhang et al. (Zhang et al., 1998; Zhang and Miller et a/., 1999) proposed an approach based on the notion of X- drop, a region that scores below -X for some fixed X>0. They consider only alignments that contain no X-drop. A more recent approach to fix the problems is the use of normalized scores for local alignments (Marzal and Vidal 1993; Arslan and 6. Egecioglu (1999); Arslan et a/., 2001). Arslan eta/. (2001) used fractional programming to obtain the optimal normalized local alignment of two sequences. They also extended their solution to find all local alignments whose normalized scores are larger than a given threshold by running their solution repeatedly after masking those alignments that are already found. In this paper we consider the problem of finding all local alignments whose normalized scores are above a given threshold. A local alignment whose normalized score is above a threshold will be called a TNL alignment. Since small regions are reported in the local alignment problem, finding all of TNL alignments can be more appropriate than just finding a single (optimal) local alignment in applications such as gene prediction (Gusfield, 1997). We present a fast algorithm for the problem of finding all the TNL alignments. Our algorithm is based on fractional programming and the banded Smith-Waterman algorithm. The fractional programming approach uses the Smith- Waterman algorithm repeatedly to obtain one local alignment. Arslan et al. (Arslan et al., 2001) use this solution repeatedly to find all the TNL alignments. Hence, Arslan et al.' s algorithm makes a double repetition of the Smith-Waterman algorithm, which makes it very slow in practice. Our algorithm first makes a careful selection of bands based on several parameters so that the selected bands can only contain TNL alignments with high probability. Then it runs the Smith-Waterman algorithm on the selected bands in such a way that the number of applications of the Smith-Waterman algorithm on each band is much smaller than that of Arslan etal.' s. Experiments show that our algorithm is about 180-330 times faster than Arslan et al.' s for the data set of human chromosome X cosmids Qc8D3 (GenBank Acc. No. AF030876 (113 kbp)) and mouse chromosome X clone BAC B22804 (GenBank Acc. No. AF121351 (123 kbp)) and it is about 40-50 times faster for the data set of human (GenBank Acc. No. L10347 (31 kbp)) and mouse (GenBank Acc. No. M65161 (32.4 kbp)) pro-alpha1 type the benefit of II collagen, while finding all TNL alignments. In addition to extremely fast speed, our algorithm can change the parameters for selecting bands to adjust the time and accuracy of our algorithm. A string or a sequence is concatenations of zero or more characters from an alphabet X. A space is denoted by J £ 2; we regard J as a character for convenience. The length of a string a is denoted by |a|. Let a, denote /th character of a string a and a(i, j) denote a substring aa+1 ■■■aj of a. When a sequence a is a substring of a sequence a, we denote it by «<a. Given two strings a=aici2---a,, and b=bib2 --bn, an alignment of a and b is a*=a2*a2*--a* and b*=hi*fe*- fo* constructed by inserting zero or more J's into a and b so that each a* maps to fo* for 1 </<I. There are three kinds of mappings in a* and b* according to the characters of and a*and b*. ■ match : a*=bi* d, • mismatch : (a* =Ab,*) and (a* bi* =/= J), • insertion or deletion (indel for short): either a* or b* is J. Note that we do not allow the case of «*=&*= J. If there exists an alignment of a and b whose number of matches, mismatches and indels are x,y and z, respectively, then we call the triplet (x, y, z) an alignment vector of a and b. We define the alignment score of an alignment a* and b* with alignment vector (x, y, z) by score(x, y, z) =x-8y ~nz that is, we assume that a match score is 1, a mismatch penalty is 8 and an indel penalty is y. We define the similarity of two sequences a and b, denoted by S(a, b), to be the highest alignment score of all possible alignments of a and b. Given two sequences a and b, a local alignment of a and b is an alignment of two strings a and p such a(a and ft{b, and the optimal local alignment of a and b is the local alignment of a and b that has the highest similarity among all local alignments of a and b. We denote the similarity of an optimal local alignment by SL(a, b). A most well-known algorithm to find an optimal local alignment was given by Smith and Waterman, which is known as the Smith-Waterman algorithm (SW algorithm for short) (Smith and Waterman, 1981; Waterman, 1995). Given two sequences a(\a\=m) and b(\b\=n), the SW algorithm computes SL(a, b) using a dynamic programming table (called the H table) of size (m+1)(n+1). Let H., for ()<i <m and 0<j<n denote SL(a(l, i), b(l,j)). Then, can be computed by the following recurrence: Among all Hu for 0<i<m and 0<j<n, the highest value is SL(a, b), and the SW algorithm takes 0(mn) time to compute it. Since the SW algorithm does not consider the lengths of local alignments, the solution of the SW algorithm may include some regions with very low similarity. Moreover, a biologically important short alignment may not be detected (Arslan and Pevzner et al., 2001; Alexandrov and Solovyev, 1998; Zhang and Miller etal., 1999). To handle these problems, the notion of normalization for similarity has been introduced (Arslan 6. Egecioglu (1999), Arslan et al., 2001; Egecioglu and Ibel, 1996; Marzal and Vidal, 1993; Alexandrov and Solovyev, 1998). Given two sequences a and b, we define the normalized similarity of a and b by, S(a, b)/(\aMb\+L), where L is a constant. From this definition we can define a normalized alignment score for local alignment: the normalized alignment score of an alignment a* and b* with alignment vector (x, y, z) is defined by nscore (x, y, z) = \+2 y +z+L (Arslan and Pevzner etal., 2001). The opti­ ma/ normalized local alignment of a and b is the local alignment of a and b that has the highest normalized alignment score. Here we consider two problems related to normalized local alignment: one is to find the optimal normalized local alignment of a and fr, the other is to find all normalized local alignments of a and b above some threshold Ts. Given two sequences a and b, let AV(a, b) be the set of all possible alignment vectors of a{a that ft{b. The problem of finding an optimal local alignment (LA problem for short) and that of finding an optimal normalized local alignment (NLA problem for short) are defined as follows: LA problem : Find (x, y, z) e AV(a, b) such that score(x, y, z) is greatest inAVfaz, b). NLA problem : Find (x, y, z) e AV(a, b) such that nscore(x, y, z) is greatest in AV(a, b). Note that we can find the score of an optimal local alignment using the SW algorithm but with a simple modification of storing the starting position of each score, we can also find the location of the optimal local alignment. For a given A, we define the parametric local alignment problem as follows: LA(A) problem : For all (x, y, z)EAV(a, b), find the maximal value ofx-J* y - //z-^(2x+2y+z+L) A parametric local alignment problem can be converted into a local alignment problem because the optimal solution of LA(A) involves solving a local alignment problem and performing some simple computations (Arslan and Pevzner etal., 2001). We can use Dinkelbach’s algorithm to solve the NLA problem. Dinkelbach developed a general algorithm that uses fractional programming (Dinkelbach, 1967). The optimal solution to NLA can be obtained via a series of optimal solutions of LA(A) for different A’s using the following result: A is the optimal solution for NLA if and only if LA(A)=0. Dinkelbach’s algorithm (Dinkelbach, 1967) (See Algorithm 1) starts with an initial value for A and repeatedly solves LA(A), Since the LA(A) problem can be converted into an LA problem, it solves the LA problem and obtains an optimal alignment vector (x, y, z). Next, it sets A= nscore(x, v, z) and repeats this process until A is not changed any more. The convergence to the optimal solution of the NLA problem is guaranteed (Arslan et al., 2001). In some biological applications, we need to find all local alignments whose normalized alignment scores are above a given threshold. Formally, we define this problem as follows: TNLA problem : For given a threshold value Ts, find all (x, y, z) e AV(a, b) such that nscore(x, y, z) >Ts. A local alignment whose normalized alignment score is above a given threshold will be called a TNL alignment. Therefore, the TNLA problem is to find all the TNL alignments. For the TNLA problem, Arslan et al. (2001) first solve the NLA problem and then mask the solution to repeatedly find the next optimal solution of the NLA problem. In general, it is highly likely that an optimal local alignment has a long part of exact matches in it. The banded Smith- Waterman algorithm uses this phenomenon and finds a solution very quickly for the LA problem with high probability (Setubal and Meidanis, 1997). Many biological systems such as Phrap (Green), STROLL (Chen and Skiena, 1997) and FASTA (Lipman and Pearson, 1998) are based on this algorithm. The banded SW algorithm consists of the following two steps: 1. Determine bands Find all the bands that can have local alignments with high similarity. Usually, a band is defined by the information of exact matches and if two or more bands overlap, they are merged into one band. 2. Run the SW algorithm Run the SW algorithm on the defined bands. The entries of Hi.j outside the bands are assumed 0 and they are not computed. We first present an algorithm for finding an optimally normalized local alignment (the NLA problem) and then present algorithms for finding all local alignments whose normalized alignment scores are above a given threshold Ts (the TNLA problem). Our algorithm for the NLA problem is used as a subroutine in our algorithms for the TNLA problem. We present an algorithm for the NLA problem of two strings a(\a\=m) and b(\b\=ri). Our algorithm is based on Dinkelbach’s algorithm which is a general algorithm to solve the NLA problem. We use the banded SW algorithm to solve the LA problem in Dinkelbach’s algorithm. Since we already described Dinkelbach’s algorithm and the banded SW algorithm, we concentrate on describing how to determine bands. We first introduce some definitions before we describe how to determine the bands. Consider the H table generated by the SW algorithm in computing SLfa, b). The fth-diagonal (or diagonal i), Q<i<m(resp. -n<f<-1), represents a diagonal that passes through (0, i) (resp. (-i, 0)) element of the H table. A band is a set of diagonals and the width of a band is the number of diagonals in the band. A band i of width k is a set of diagonals from i through i+k- 1. We say that an exact match (i,j) occurs in diagonal k if a(i,j)= b(i+k,j+k). Exact match (i,j) is maximal if a-i^bi+k-i and ai+i=Abj+M. The length of exact match (i, j) is j-i+1. We show how to determine the bands with given parameters 7), w, and 7k The parameter T is a threshold on the length of an exact match, w is the width of a band, and Tb is a threshold on the weight of a band. Determining the bands is composed of the following four steps: 1. Find all maximal exact matches longer than T in each diagonal. We first generate all suffixes of a and b and sort them in lexicographic order. Once all the suffixes of a and b are sorted, it is easy to find all maximal exact matches longer than T. 2. Compute the weight of each diagonal, where the weight of a diagonal is the sum of the lengths of all maximal exact matches longer than T in the diagonal. 3. Select every band b of width w whose weight is above Tb, where the weight of a band is the sum of the weights of all diagonals in the band. 4. Merge two bands Bi={di, dj} and B 2 ={dk, — di} into a single band B={min(d, dk), •••, mwcfd, di)}, if they overlap (i.e., i<k<j or i<l<j). Repeat merging until no two bands overlap. (Fig. 1) Parameters T, w, and Tb should be chosen appropriately so that an optimal normalized local alignment does not lie out of the bands. In Section 4, we suggest appropriate values for Ti, w, and Tb through experiments. Our algorithm for the NLA problem first determines the bands as explained above and then runs the SW algorithm on the selected bands. The speedup achieved by our algorithm over Arslan et al.' s is due to the following two improvements. ■ We uses the banded SW algorithm to solve the LA problem of Dinkelbach’s algorithm while Arslan et al. used the SW algorithm. Thus, the speedup achieved by the banded SW algorithm over the SW algorithm is reflected in our algorithm. We perform the step of determining bands in the banded SW algorithm only once while we perform the banded SW algorithm several times in Dinkelbach’s algorithm. Since the bands are the sets of diagonals where long exact matches occur, the bands in all repetitions of the banded SW algorithm are the same. Hence, we need to determine the bands only once. We now present algorithms for the TNLA problem of two strings a and b. We first give a simple algorithm that finds all normalized local alignments whose alignment scores are above a given threshold L (the TNL alignments) in such a way that a TNL alignment of larger alignment score is found earlier than a TNL alignment of smaller alignment score. Arslan etal. (Arslan etal., 2001) also found the TNL alignments in the same decreasing order. Then, we improve the simple algorithm by finding the TNL alignments in a somewhat different order. The simple algorithm consists of two steps. Step 1 is performed once and step 2 is repeated until all TNL alignments are found. Step 1. Find all maximal exact matches of a and b longer than Ti and determine the bands as we did in Section 3.1. Step 2. Find an optimal normalized local alignment of a and b using the algorithm presented in Section 3.1. If the alignment score of the found alignment is above T, mask the found alignment and repeat step 2. Otherwise, terminate. The simple algorithm finds the normalized local alignments in decreasing order of alignment scores. Hence, it is guaranteed that every TNL alignment has been found when the simple algorithm terminates. We consider the speedup achieved by the simple algorithm over Arslan et a/.’ s algorithm. Let p denote the number of the TNL alignments. Step 2 iterates just p times. Since the time for masking an alignment is negligible and our algorithm performs step 1 only once, our simple algorithm for the TNLA problem takes about p times our algorithm for the NLA problem. Also, Arslan et a/.’ s algorithm for the TNLA problem takes p times their algorithm for the NLA problem. Therefore, the speedup achieved by our simple algorithm over Arslan et al’ s algorithm for the TNLA problem is the speedup achieved by our algorithm over Arslan et a/.’s algorithm for the NLA problem. As long as all the TNL alignments are guaranteed to be found, the order of the alignments that are found may not be important. If we find the TNL alignments in a different order, we can find them more efficiently. The improved algorithm is as follows. Step 1. Find all maximal exact matches of a and b longer than T/ and determine the bands as in Section 3.1. Let Bi, B2, “,Bk denote the bands. Step 2. For each band B> for 1 <i<k, perform the following. Find an optimal normalized local alignment within Bi by running the SW algorithm only in Bi. If the normalized alignment score of the found alignment is above L, mask the found alignment and repeat. Otherwise, we are done with Bi. It is easy to see that all the TNL alignments have been found when the improved algorithm terminates. We consider the speedup achieved by the improved algorithm over the simple algorithm. Let C, l<i<k, denote the number of the TNL alignments included in band B. and W the width of B,. In the improved algorithm, band B, accessed when we find each of C TNL alignments while the simple algorithm, all bands Bi, B2, •••, B k are accessed when we find each of G+C2+—+G TNL alignments. Note that the improved algorithm accesses B, at least once even if there is no TNL alignment in it. Hence, the speedup achieved by the improved algorithm over the simple algorithm is Remark: We can consider affine gaps in both problems of NLA and TNLA. The affine gap penalty for a gap is defined as 7+/*k, where 7 is a gap open penalty, p is an indel penalty and k is the length of a gap. We can include the affine gap penalty in our algorithms for NLA and TNLA by a simple modification of the recurrence of the SW algorithm (Gotoh 1982; Waterman 1995). We implemented and compared two algorithms for the NLA problem: Arslan et al.' s and ours, and three algorithms for the TNLA problem: Arslan et al’.s algorithm, our simple algorithm and our improved algorithm. In all our implementations, we considered the affine gaps. We have chosen two data sets to test the efficiency of the algorithms: (i) human chromosome X cosmids Qc8D3 (GenBank Acc. No. AF030876 (113 kbp)) and mouse chromosome X clone BAC B22804 (GenBank Acc. No. AF121351 (123 kbp)) and (ii) human (GenBank Acc. No. L10347 (31 kbp)) and mouse (GenBank Acc. No. M65161 (32.4 kbp)) pro-alpha1 type II collagen. Since the repeats which are biologically uninteresting regions may have high normalized scores, we masked the repeats by RepeatMasker (http:// repeatmasker, genome.washington. edu/) before running all the algorithms. We implemented all the algorithms in C++ language. The programs were run on Pentium III 866MHzx2 work­ station with 768MB main memory. The parameter values used in our algorithms are as follows: L=2000 (a constant in defining normalized scores), £=1 (mismatch penalty), 7 =6 (gap open penalty), /z=0.2 (indel penalty), and T=0.035 (threshold of normalized score). Our algorithm finds the same optimal normalized local alignment as Arslan et al.’s algorithm does even when w=100 for both data sets. As shown in Table 1, our algorithm is about 20 times faster than Arslan et al.' s for both data sets when w> 150 and Ti< 150. We have experimented with various values of the three parameters Ti, T> and w, and here we show nine cases for three values of Ti and three values of w and we have chosen an appropriate value of TI for each T. Note that we decrease T b as Ti increases. Otherwise, good bands may be discarded since longer exact matches are fewer than shorter exact matches. Table 2 and Table 3 show execution times of the three algorithms for the TNLA problem. We first analyze the ratios of execution times of the three algorithms. Assume that the widths of all bands are the same and the TNL alignments are equally distributed in the bands which include TNL alignments. Let k be the number of the bands that our algorithms determine, and k’ be the number of the bands which include TNL alignments. Also, let t be the average number of TNL alignments in a band. Then, the equation (1) can be rewritten as follows: kk't k't+k-k' Hence, our improved algorithm is faster than Arslan et al’s by the following factor: where 5 is the speedup of our simple algorithm over Arslan et al.’ s. Consider the experimental results for the first data set in Table 2. When w=150, 7/=15 and T>=50, we have 5=80 and k=5, k’=5, t=\/3 (k, k’ are not shown in Table 2 since the number of TNL alignments is 15 as shown in Table 4). By the above formula, our improved algorithm should be approximately 400 times faster than Arslan et a/.’s and the result (631603/1895=330) is similar to our expectation. (Note that the time of the improved algorithm can vary depending on the value of 71 in each case.) For another example, when w=200, T<=12 and 71=80, we get 5=18, £=24, £ =5 and r=1/3. Thus, our improved algorithm should be about 200 times faster than Arslan etal’ s and the result (631603/3492=180) shows an approximately same ratio. For the second data set, our algorithm is about 40 times faster than Arslan et al’ s when w=150, 77=15 and 71=50, and about 50 times faster when w=200, 77=12 and 71=80. Table 4 shows the accuracy of our algorithms. For example, when w=100, 77=20 and 71=30 in Table 4, our algorithms find 12 out of 15 TNL alignments. But as w increases and T decreases, our algorithms are getting more accurate. Note that our algorithms find all TNL alignments when w> 150 and T< 15. In addition to extremely fast speeds, another advantage of our algorithms is that we can change the parameters w, Ti and 71 to adjust the time and accuracy of our algorithms.