1
2 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
3
4 To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,dflag,&SI,&SJ,&EI,&EJ,MW).
5 The parameters are explained as follows:
6 A, B : two sequences to be aligned
7 M : the length of sequence A
8 N : the length of sequence B
9 L : lower bound of the band
10 U : upper bound of the band
11 W : scoring table for matches and mismatches
12 G : gap-opening penalty
13 H : gap-extension penalty
14 dflag : 0 - no display or backward pass
15 *SI : starting position of sequence A in the optimal local alignment
16 *SJ : starting position of sequence B in the optimal local alignment
17 *EI : ending position of sequence A in the optimal local alignment
18 *EJ : ending position of sequence B in the optimal local alignment
19 MW : maximum window size
20 */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #ifdef BIGMEM
26 #define MININT -9999999
27 #else
28 #define MININT -32000
29 #endif
30
31 static int *CC=NULL, *DD;
32
33 #define max(x,y) ((x) >= (y) ? (x) : (y))
34 #define min(x,y) ((x) <= (y) ? (x) : (y))
35
LOCAL_ALIGN(A,B,M,N,low,up,W,G,H,dflag,psi,psj,pei,pej,MW)36 int LOCAL_ALIGN(A,B,M,N,low,up,W,G,H,dflag,psi,psj,pei,pej,MW)
37 char A[],B[]; int M,N,low,up; int **W, G,H;
38 int dflag;
39 int *psi, *psj, *pei, *pej;
40 int MW;
41 {
42 int band;
43 char *ckalloc();
44 int i, j, si, ei;
45 int c, d, e, t, m;
46 int leftd, rightd;
47 int best_score, starti, startj, endi, endj;
48 int *wa, curd;
49 int ib;
50 char flag;
51
52 m = G+H;
53 low = max(-M, low);
54 up = min(N, up);
55
56 if (N <= 0) {
57 *psi = *psj = *pei = *pej;
58 return 0;
59 }
60 if (M <= 0) {
61 *psi = *psj = *pei = *pej;
62 return 0;
63 }
64 band = up-low+1;
65 if (band < 1) {
66 printf("low > up is unacceptable!\n");
67 exit(1);
68 }
69 j = (MW+2+2) * sizeof(int);
70 if (CC==NULL) {
71 CC = (int *) ckalloc(j);
72 DD = (int *) ckalloc(j);
73 }
74
75 if (low > 0) leftd = 1;
76 else if (up < 0) leftd = band;
77 else leftd = 1-low;
78 rightd = band;
79 si = max(0,-up);
80 ei = min(M,N-low);
81 CC[leftd] = 0;
82 for (j = leftd+1; j <= rightd; j++) {
83 CC[j] = 0;
84 DD[j] = -G;
85 }
86 CC[rightd+1] = MININT;
87 DD[rightd+1] = MININT;
88 best_score = 0;
89 endi = si;
90 endj = si+low;
91 CC[leftd-1] = MININT;
92 DD[leftd] = -G;
93 for (i = si+1; i <= ei; i++) {
94 if (i > N-up) rightd--;
95 if (leftd > 1) leftd--;
96 wa = W[A[i]];
97 if ((c = CC[leftd+1]-m) > (d = DD[leftd+1]-H)) d = c;
98 if ((ib = leftd+low-1+i ) > 0) c = CC[leftd]+wa[B[ib]];
99 /*
100 if (ib > N) fprintf(stderr,"B[%d] out of range %d\n",ib,N);
101 */
102 if (d > c) c = d;
103 if (c < 0) c = 0;
104 e = c-G;
105 DD[leftd] = d;
106 CC[leftd] = c;
107 if (c > best_score) {
108 best_score = c;
109 endi = i;
110 endj = ib;
111 }
112 for (curd=leftd+1; curd <= rightd; curd++) {
113 if ((c = c-m) > (e = e-H)) e = c;
114 if ((c = CC[curd+1]-m) > (d = DD[curd+1]-H)) d = c;
115 /*
116 if ((ib=curd+low-1+i) <= 0 || ib > N)
117 fprintf(stderr,"B[%d]:%d\n",ib,B[ib]);
118 */
119 c = CC[curd] + wa[B[curd+low-1+i]];
120 if (e > c) c = e;
121 if (d > c) c = d;
122 if (c < 0) c = 0;
123 CC[curd] = c;
124 DD[curd] = d;
125 if (c > best_score) {
126 best_score = c;
127 endi = i;
128 endj = curd+low-1+i;
129 }
130 }
131 }
132
133 if (!dflag) {
134 *pei = endi;
135 *pej = endj;
136 return best_score;
137 }
138
139 leftd = max(1,-endi-low+1);
140 rightd = band-(up-(endj-endi));
141 CC[rightd] = 0;
142 t = -G;
143 for (j = rightd-1; j >= leftd; j--) {
144 CC[j] = t = t-H;
145 DD[j] = t-G;
146 }
147 for (j = rightd+1; j <= band; ++j) CC[j] = MININT;
148 CC[leftd-1] = DD[leftd-1] = MININT;
149 DD[rightd] = -G;
150 flag = 0;
151 for (i = endi; i >= 1; i--) {
152 if (i+low <= 0) leftd++;
153 if (rightd < band) rightd++;
154 wa = W[A[i]];
155 if ((c = CC[rightd-1]-m) > (d = DD[rightd-1]-H)) d = c;
156 if ((ib = rightd+low-1+i) <= N) c = CC[rightd]+wa[B[ib]];
157
158 /*
159 if (ib <= 0) fprintf(stderr,"rB[%d] <1\n",ib);
160 */
161 if (d > c) c = d;
162 e = c-G;
163 DD[rightd] = d;
164 CC[rightd] = c;
165 if (c == best_score) {
166 starti = i;
167 startj = ib;
168 flag = 1;
169 break;
170 }
171 for (curd=rightd-1; curd >= leftd; curd--) {
172 if ((c = c-m) > (e = e-H)) e = c;
173 if ((c = CC[curd-1]-m) > (d = DD[curd-1]-H)) d = c;
174
175 /*
176 if ((ib=curd+low-1+i) <= 0 || ib > N)
177 fprintf(stderr,"i: %d, B[%d]:%d\n",i,ib,B[ib]);
178 */
179 c = CC[curd] + wa[B[curd+low-1+i]];
180 if (e > c) c = e;
181 if (d > c) c = d;
182 CC[curd] = c;
183 DD[curd] = d;
184 if (c == best_score) {
185 starti = i;
186 startj = curd+low-1+i;
187 flag = 1;
188 break;
189 }
190 }
191 if (flag == 1) break;
192 }
193
194 if (starti < 0 || starti > M || startj < 0 || startj > N) {
195 printf("starti=%d, startj=%d\n",starti,startj);
196 *psi = *psj = *pei = *pej;
197 exit(1);
198 }
199 *psi = starti;
200 *psj = startj;
201 *pei = endi;
202 *pej = endj;
203 return best_score;
204 }
205
206 #ifdef LFASTA
207
LLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,dflag,psi,psj,pei,pej,MW)208 int LLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,dflag,psi,psj,pei,pej,MW)
209 char A[],B[];
210 int M,N,low,up;
211 int **W,G,H;
212 int dflag;
213 int *psi, *psj, *pei, *pej;
214 int MW;
215 {
216 int band;
217 char *ckalloc();
218 int i, j, si, ei;
219 int c, d, e, t, m;
220 int lc; /* row local score */
221 int leftd, rightd;
222 int best_score, starti, startj, endi, endj;
223 int *wa, curd, ib;
224 int flag;
225
226 m = G+H;
227 low = max(-M, low);
228 up = min(N, up);
229 if (N <= 0) {
230 *psi = *psj = *pei = *pej;
231 return 0;
232 }
233 if (M <= 0) {
234 *psi = *psj = *pei = *pej;
235 return 0;
236 }
237 band = up-low+1;
238 if (band < 1) {
239 printf("low > up is unacceptable!\n");
240 exit(1);
241 }
242
243 j = (MW+1+2) * sizeof(int);
244 if (CC==NULL) {
245 CC = (int *) ckalloc(j);
246 DD = (int *) ckalloc(j);
247 }
248
249 if (low > 0) leftd = 1;
250 else if (up < 0) leftd = band;
251 else leftd = 1-low;
252 rightd = band;
253 si = max(0,-up);
254 ei = min(M,N-low);
255 CC[leftd] = 0;
256 for (j = leftd+1; j <= rightd; j++) {
257 CC[j] = 0;
258 DD[j] = -G;
259 }
260 CC[rightd+1] = MININT;
261 DD[rightd+1] = MININT;
262 best_score = 0;
263 endi = si;
264 endj = si+low;
265 CC[leftd-1] = MININT;
266 DD[leftd] = -G;
267 for (i = si+1; i <= ei; i++) {
268 if (i > N-up) rightd--;
269 if (leftd > 1) leftd--;
270 e = MININT;
271 wa = W[A[i]];
272 if ((c = CC[leftd+1]-m) > (d = DD[leftd+1]-H)) d = c;
273 if ((ib = leftd+low-1+i) > 0) c = CC[leftd]+wa[B[ib]];
274 if (d > c) c = d;
275 if (c < 0) c = 0;
276 e = c-G;
277 DD[leftd] = d;
278 CC[leftd] = c;
279 if (c > best_score) {
280 best_score = c;
281 endi = i;
282 endj = ib;
283 }
284 for (curd=leftd+1; curd <= rightd; curd++) {
285 if ((c = c-m) > (e = e-H)) e = c;
286 if ((c = CC[curd+1]-m) > (d = DD[curd+1]-H)) d = c;
287 c = CC[curd] + wa[B[curd+low-1+i]];
288 if (e > c) c = e;
289 if (d > c) c = d;
290 if (c < 0) c = 0;
291 if (c > lc) lc = c;
292 CC[curd] = c;
293 DD[curd] = d;
294 if (c > best_score) {
295 best_score = c;
296 endi = i;
297 endj = curd+low-1+i;
298 }
299 }
300 if (lc <= 0) break;
301 }
302
303 if (!dflag) {
304 /* free(CC); free(DD); */
305 *pei = endi;
306 *pej = endj;
307 return best_score;
308 }
309 }
310
RLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,psi,psj,pei,pej,MW)311 int RLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,psi,psj,pei,pej,MW)
312 char A[],B[];
313 int M,N,low,up;
314 int **W,G,H;
315 int *psi, *psj, *pei, *pej;
316 int MW;
317 {
318 int band;
319 char *ckalloc();
320 int i, j, si, ei;
321 int c, d, e, t, m;
322 int lc; /* row local score */
323 int leftd, rightd;
324 int best_score, starti, startj;
325 int *wa, curd, ib;
326
327 m = G+H;
328 low = max(-M, low);
329 up = min(N, up);
330 if (N <= 0) {
331 *psi = *psj = *pei = *pej;
332 return 0;
333 }
334 if (M <= 0) {
335 *psi = *psj = *pei = *pej;
336 return 0;
337 }
338 band = up-low+1;
339 if (band < 1) {
340 printf("low > up is unacceptable!\n");
341 exit(1);
342 }
343 j = (MW+1+2) * sizeof(int);
344
345 if (CC==NULL) {
346 CC = (int *) ckalloc(j);
347 DD = (int *) ckalloc(j);
348 }
349
350 leftd = max(1,-M-low+1);
351 rightd = band-(up-(N-M));
352 CC[rightd] = 0;
353 t = -G;
354 for (j = rightd-1; j >= leftd; j--) {
355 CC[j] = t = t-H;
356 DD[j] = t-G;
357 }
358 CC[leftd-1] = DD[leftd-1] = MININT;
359 CC[rightd+1] = MININT;
360 DD[rightd] = -G;
361 best_score = 0;
362 starti = M;
363 startj = N;
364 for (i = M; i >= 1; i--) {
365 if (i+low <= 0) leftd++;
366 if (rightd < band) rightd++;
367 e = MININT;
368 wa = W[A[i]];
369 if ((c = CC[rightd-1]-m) > (d = DD[rightd-1]-H)) d = c;
370 if ((ib = rightd+low-1+i) <= N) c = CC[rightd]+wa[B[ib]];
371 if (d > c) c = d;
372 e = c-G;
373 DD[rightd] = d;
374 CC[rightd] = c;
375 if (c == best_score) {
376 starti = i;
377 startj = ib;
378 break;
379 }
380 lc = 0;
381 for (curd=rightd-1; curd >= leftd; curd--) {
382 if ((c = c-m) > (e = e-H)) e = c;
383 if ((c = CC[curd-1]-m) > (d = DD[curd-1]-H)) d = c;
384 c = CC[curd] + wa[B[curd+low-1+i]];
385 if (e > c) c = e;
386 if (d > c) c = d;
387 if (c < 0) c = 0;
388 if (lc < c) lc = c;
389 CC[curd] = c;
390 DD[curd] = d;
391 if (c > best_score) {
392 best_score = c;
393 starti = i;
394 startj = curd+low-1+i;
395 }
396 }
397 if (lc <= 0) break;
398 }
399
400 /* free(CC);
401 free(DD);
402 */
403 if (starti < 0 || starti > M || startj < 0 || startj > N) {
404 printf("starti=%d, startj=%d\n",starti,startj);
405 *psi = *psj = *pei = *pej;
406 exit(1);
407 }
408 *psi = starti;
409 *psj = startj;
410 return best_score;
411 }
412
ckalloc(amount)413 char *ckalloc(amount)
414 int amount;
415 {
416 char *p;
417
418 if ((p = malloc( (unsigned) amount)) == NULL)
419 fatal("Ran out of memory.");
420 return(p);
421 }
422
423 /* fatal - print message and die */
fatal(msg)424 fatal(msg)
425 char *msg;
426 {
427 fprintf(stderr, "%s\n", msg);
428 exit(1);
429 }
430
431 /* dummy for zzlgmata.c */
432
ALIGN()433 int ALIGN() {}
434
435 #endif
436
437