1
2 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
3
4 To invoke, call FLOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,MW).
5 The parameters are:
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 MW : maximum window size
15 */
16
17 #include <stdio.h>
18 #include <stdlib.h>
19
20 #ifdef BIGMEM
21 #define MININT -9999999
22 #else
23 #define MININT -32000
24 #endif
25
26 static struct swstr {
27 int CC;
28 int DD;
29 } *ss;
30
31 #ifndef max
32 #define max(x,y) ((x) >= (y) ? (x) : (y))
33 #define min(x,y) ((x) <= (y) ? (x) : (y))
34 #endif
35
FLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,MW)36 int FLOCAL_ALIGN(A,B,M,N,low,up,W,G,H,MW)
37 char A[],B[];
38 int M,N,low,up;
39 int **W,G,H;
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;
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) return 0;
57
58 if (M <= 0) return 0;
59
60 band = up-low+1;
61 if (band < 1) {
62 fprintf(stderr,"low > up (%d > %d)\n",low,up);
63 return 0;
64 }
65
66 j = (MW+2+2) * sizeof(struct swstr);
67 if (ss==NULL) ss = (struct swstr *) ckalloc(j);
68
69 if (low > 0) leftd = 1;
70 else if (up < 0) leftd = band;
71 else leftd = 1-low;
72
73 rightd = band;
74 si = max(0,-up);
75 ei = min(M,N-low);
76 ss[leftd].CC = 0;
77
78 for (j = leftd+1; j <= rightd; j++) {
79 ss[j].CC = 0;
80 ss[j].DD = -G;
81 }
82 ss[rightd+1].CC = MININT;
83 ss[rightd+1].DD = MININT;
84
85 best_score = 0;
86 ss[leftd-1].CC = MININT;
87 ss[leftd].DD = -G;
88
89 for (i = si+1; i <= ei; i++) {
90 if (i > N-up) rightd--;
91 if (leftd > 1) leftd--;
92 wa = W[A[i]];
93 if ((c = ss[leftd+1].CC-m) > (d = ss[leftd+1].DD-H)) d = c;
94 if ((ib = leftd+low-1+i ) > 0) c = ss[leftd].CC+wa[B[ib]];
95
96 if (d > c) c = d;
97 if (c < 0) c = 0;
98 e = c-G;
99 ss[leftd].DD = d;
100 ss[leftd].CC = c;
101 if (c > best_score) best_score = c;
102
103 for (curd=leftd+1; curd <= rightd; curd++) {
104 if ((c = c-m) > (e = e-H)) e = c;
105 if ((c = ss[curd+1].CC-m) > (d = ss[curd+1].DD-H)) d = c;
106 c = ss[curd].CC + wa[B[curd+low-1+i]];
107 if (e > c) c = e;
108 if (d > c) c = d;
109 if (c < 0) c = 0;
110 ss[curd].CC = c;
111 ss[curd].DD = d;
112 if (c > best_score) best_score = c;
113 }
114 }
115 return best_score;
116 }
117
118