← Back to branch summary

~ubuntu-branches/ubuntu/warty/ncbi-tools6/warty

« back to all changes in this revision

Viewing changes to tools/bandalg3.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2002-04-04 22:13:09 UTC
  • Revision ID: james.westby@ubuntu.com-20020404221309-vfze028rfnlrldct
Tags: upstream-6.1.20011220a
Import upstream version 6.1.20011220a

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/bandalg3.c
 
1
/* ===========================================================================
 
2
 *
 
3
 *                            PUBLIC DOMAIN NOTICE
 
4
 *               National Center for Biotechnology Information
 
5
 *
 
6
 *  This software/database is a "United States Government Work" under the
 
7
 *  terms of the United States Copyright Act.  It was written as part of
 
8
 *  the author's official duties as a United States Government employee and
 
9
 *  thus cannot be copyrighted.  This software/database is freely available
 
10
 *  to the public for use. The National Library of Medicine and the U.S.
 
11
 *  Government have not placed any restriction on its use or reproduction.
 
12
 *
 
13
 *  Although all reasonable efforts have been taken to ensure the accuracy
 
14
 *  and reliability of the software and data, the NLM and the U.S.
 
15
 *  Government do not and cannot warrant the performance or results that
 
16
 *  may be obtained by using this software or data. The NLM and the U.S.
 
17
 *  Government disclaim all warranties, express or implied, including
 
18
 *  warranties of performance, merchantability or fitness for any particular
 
19
 *  purpose.
 
20
 *
 
21
 *  Please cite the author in any work or product based on this material.
 
22
 *
 
23
 * ==========================================================================*/
 
24
 
 
25
/****************************************************************************
 
26
 
 
27
  File name: bandalg3.c
 
28
  
 
29
  Author: Gennadiy Savchuk, Jinqhui Zhang, Tom Madden
 
30
  
 
31
  Contents: Functions to perform global banded alignment.
 
32
 
 
33
  ***************************************************************************/
 
34
 
 
35
#include <bandalgn.h>
 
36
 
 
37
static Int4 g_band3_align(Uint1Ptr A, Uint1Ptr B,
 
38
                          Int4 M, Int4 N,
 
39
                          Int4 low, Int4 up, data_t *data);
 
40
 
 
41
static Int4 g_band3_CHECK_SCORE(Uint1Ptr A, Uint1Ptr B,
 
42
                                Int4 M, Int4 N,
 
43
                                Int4 PNTR S, data_t *data);
 
44
 
 
45
Int4 LIBCALL gband_linear_qgap(Uint1Ptr A, Uint1Ptr B,
 
46
                               Int4 M, Int4 N,
 
47
                               Int4Ptr PNTR matrix,
 
48
                               PSUGapOptionsPtr option,
 
49
                               Int4Ptr S, Int4Ptr Slen)
 
50
{ 
 
51
  data_t data;
 
52
  Int4 c, i, j;
 
53
  Int4 low, up;
 
54
  Int4 band;
 
55
  Int4 score;
 
56
 
 
57
  /* Setup global parameters */
 
58
  data.g = option->gopen;
 
59
  data.zzh = option->gext;
 
60
  data.w = matrix;
 
61
  data.m = data.g + data.zzh;
 
62
 
 
63
  data.sapp = S;
 
64
  data.last = 0;
 
65
  *Slen = 0;
 
66
 
 
67
  low = option->start_diag;
 
68
  band = option->width;
 
69
  up = band + low - 1;
 
70
 
 
71
  low = MIN(MAX(-M, low),MIN(N-M,0));
 
72
  up = MAX(MIN(N, up),MAX(N-M,0));
 
73
 
 
74
  if(up < 0 || low > 0) {
 
75
    ErrPostEx(SEV_WARNING, 0, 0,
 
76
              "The band does not include (0,0) grid point");
 
77
    return 0;
 
78
  } 
 
79
  if(up+M < N || low+M > N) {
 
80
    ErrPostEx(SEV_WARNING, 0, 0,
 
81
              "The band does not include lower corner");
 
82
    return 0;
 
83
  }
 
84
  
 
85
  if(N <= 0) { 
 
86
    if(M > 0) DEL_(M);
 
87
    return -gap_(M);
 
88
  }
 
89
  if(M <= 0) {
 
90
    INS_(N);
 
91
    return -gap_(N);
 
92
  }
 
93
  if((band = up - low + 1) <= 1) {
 
94
    c = 0;
 
95
    for(i = 1; i <= M; i++) {
 
96
      REP_;
 
97
      c += data.w[A[i]][B[i]];
 
98
    }
 
99
    return c;
 
100
  }
 
101
 
 
102
  j = (band+2) * sizeof(Int4);
 
103
  data.leggA = data.leggB = data.reggA = data.reggB = data.leghA =
 
104
    data.leghB = data.reghA = data.reghB = 0;
 
105
 
 
106
  data.CD = (dp_ptr)MemGet(sizeof(dp_node) * (band + 2), MGET_ERRPOST);
 
107
 
 
108
  /* if((c = g_band3_align(A,B,M,N,low,up,0,0)) != (score = g_band3_CHECK_SCORE(A,B,M,N,S))) */
 
109
  c = g_band3_align(A, B, M, N, low, up, &data);
 
110
  score = g_band3_CHECK_SCORE(A, B, M, N, S, &data);
 
111
  if(c != score) { 
 
112
    ErrPostEx(SEV_WARNING, 0, 0, "Check_Score = %ld align_score = %ld ", 
 
113
              (long) score, (long) c);
 
114
        return 0;
 
115
  }
 
116
  MemFree(data.CD);
 
117
 
 
118
  *Slen = data.sapp - S;
 
119
 
 
120
  return c;
 
121
}
 
122
 
 
123
/* g_band3_CHECK_SCORE - return the score of the alignment stored in S */
 
124
 
 
125
static Int4 g_band3_CHECK_SCORE(Uint1Ptr A, Uint1Ptr B,
 
126
                                Int4 M, Int4 N,
 
127
                                Int4 PNTR S, data_t *data)
 
128
{ 
 
129
  register Int4 i, j, op;
 
130
  Int4 score;
 
131
 
 
132
  score = i = j = op = 0;
 
133
  if (*S < 0) {
 
134
    score -= data->leggA - *S * data->leghA;
 
135
    i = i - *S++;
 
136
  } else if ((*S) > 0) {
 
137
    score -= data->leggB + *S * data->leghB;
 
138
    j = j + *S++;
 
139
  }
 
140
  while (i < M || j < N) {
 
141
    op = *S++;
 
142
    if(op == 0) 
 
143
      score = data->w[A[++i]][B[++j]] + score;
 
144
    else if(op > 0) {
 
145
      score = score - (data->g + op * data->zzh);
 
146
      j = j + op;
 
147
    }
 
148
    else {
 
149
      score = score - (data->g - op * data->zzh);
 
150
      i = i - op;
 
151
    }
 
152
  }
 
153
  if(op < 0)
 
154
    score += (data->g - op * data->zzh) - (data->reggA - data->reghA * op);
 
155
  else if(op > 0)
 
156
    score += (data->g + op * data->zzh) - (data->reggB + data->reghB * op);
 
157
 
 
158
 
 
159
  return(score);
 
160
}
 
161
 
 
162
 
 
163
/* g_band3_align(A, B, M, N, up, low, tb, te, data) returns the cost
 
164
   of an optimum conversion between
 
165
   A[1..M] and B[1..N] and appends such a conversion to the current script.
 
166
   tb(te)= 1  no gap-open penalty if the conversion begins(ends) with a delete.
 
167
   tb(te)= 2  no gap-open penalty if the conversion begins(ends) with an insert.
 
168
*/
 
169
static Int4 g_band3_align(Uint1Ptr A, Uint1Ptr B,
 
170
                          Int4 M, Int4 N,
 
171
                          Int4 low, Int4 up, data_t *data)
 
172
{
 
173
  Int4 k, v;
 
174
  Int4 band, j;
 
175
  Int4 leftd, rightd;   /* for CC, DD, CP and DP */
 
176
  register Int4 curd;   /* current index for CC, DD CP and DP */
 
177
  register Int4 i;
 
178
  register Int4 c, d, e, x;
 
179
  register dp_ptr ap;
 
180
  Int4 t;
 
181
  Int4Ptr wa;
 
182
  Int1Ptr PNTR state, st, tmp;
 
183
  Int4 ib, best=MININT, X, Y;
 
184
 
 
185
 
 
186
  /* Boundary cases: M <= 0 , N <= 0, or up-low <= 0 */
 
187
  band = up - low + 1;
 
188
  state = (Int1Ptr PNTR) MemGet(sizeof(Int1Ptr)*(M+1), MGET_ERRPOST);
 
189
  state[0] = (Int1Ptr) MemGet((M+1)*(band+2), MGET_ERRPOST);
 
190
  for (i = 1; i <= M; i++) state[i] = state[i-1]+band+2;
 
191
 
 
192
  /* Initialization */
 
193
  leftd = 1-low;
 
194
  rightd = up-low+1;
 
195
 
 
196
  data->CD[leftd].CC = 0; state[0][leftd] = -1;
 
197
 
 
198
  t = -data->leggB;
 
199
  for(j = leftd + 1; j <= rightd; j++) {
 
200
    data->CD[j].CC = t = t - data->leghB;
 
201
    data->CD[j-1].DD = t - data->m;
 
202
    state[0][j] = 1;
 
203
  }
 
204
  data->CD[rightd+1].CC = MININT;
 
205
  data->CD[rightd].DD = MININT;
 
206
  data->CD[leftd-1].DD = -data->leggA;
 
207
  data->CD[leftd-1].CC = MININT;
 
208
  for (i = 1; i <= M; i++) {
 
209
    if (i > N-up) rightd--;
 
210
    if (leftd > 1) leftd--;
 
211
    wa = data->w[A[i]];
 
212
    d = data->CD[leftd].DD;
 
213
    k = 0;
 
214
    if ((ib = leftd+low-1+i) > 0) c = data->CD[leftd].CC+wa[B[ib]];
 
215
    if (d > c || ib <= 0) {
 
216
      c = d;
 
217
      k = 2;
 
218
    }
 
219
    e = c - data->m;
 
220
    if(ib <= 0) {
 
221
      data->CD[leftd - 1].DD = d - data->leghA;
 
222
      k += 20;
 
223
    }
 
224
    st = &state[i][leftd];
 
225
    *st++ = (Int1) k;
 
226
    data->CD[leftd].CC = c;
 
227
    for(curd = leftd + 1, ap = &data->CD[curd]; curd <= rightd; curd++) {
 
228
      c = ap->CC + wa[B[curd + low - 1 + i]];
 
229
      if((d = ap->DD) > c) {
 
230
        if(d > e) {
 
231
          ap->CC = d;
 
232
          *st++ =32;
 
233
        }
 
234
        else {
 
235
          ap->CC = e;
 
236
          *st++=31;
 
237
        }
 
238
        e -= data->zzh;
 
239
        (ap++ - 1)->DD = d - data->zzh;
 
240
      }
 
241
      else if (e > c) {                
 
242
        ap->CC = e;
 
243
        e -= data->zzh;
 
244
        (ap++ - 1)->DD = d - data->zzh;
 
245
        *st++ = 31;
 
246
      }
 
247
      else {
 
248
        ap->CC = c;
 
249
        if((c -= data->m) > (e -= data->zzh)) {
 
250
          if(c > (d -= data->zzh)) {
 
251
            (ap++ - 1)->DD = e = c;
 
252
            *st++ = 0;
 
253
          }
 
254
          else {
 
255
            e = c;
 
256
            (ap++ - 1)->DD = d;
 
257
            *st++ = 20;
 
258
          } 
 
259
        }
 
260
        else {
 
261
          if(c > (d -= data->zzh)) {
 
262
            (ap++ - 1)->DD = c;
 
263
            *st++ = 10;
 
264
          }
 
265
          else {
 
266
            (ap++ - 1)->DD = d;
 
267
            *st++ = 30;
 
268
          }
 
269
        }
 
270
      }
 
271
    }
 
272
    if(i > N-up &&
 
273
       best < (j = (ap - 1)->CC - data->reggA - (N - i) * data->reghA)) {
 
274
      best = j; X = i; Y = rightd;
 
275
    }
 
276
  }
 
277
  for(ap = &data->CD[leftd]; ap <= &data->CD[rightd]; ap++) {
 
278
    if((j = ap->CC - data->reggB - data->reghB *
 
279
        (x = (&data->CD[rightd] - ap)))
 
280
       > best) {
 
281
      X = M;
 
282
      best = j;
 
283
      Y = rightd - x;
 
284
    }
 
285
  }
 
286
  if (data->CD[rightd].CC > best) {X=M; Y= N; best = data->CD[rightd].CC;}
 
287
  v= best;
 
288
  tmp = MemGet(M+N, MGET_ERRPOST);
 
289
  for (i = X, j = Y, e=0, c = 0; i>=0; i--, c++) {
 
290
    k  = (t=state[i][j]) %10;
 
291
    if (t == -1) break;
 
292
    if (e == 1 && (t/10)%2 == 1) k = 1;
 
293
    if (e == 2 && (t/20)== 1) k = 2;
 
294
    if (k == 1) { j--; i++;}
 
295
    else if (k == 2) j++;
 
296
    e = k;
 
297
    tmp[c] = (Int1) e;
 
298
  }
 
299
  for (i = c-1; i >= 0; i--) 
 
300
    switch(tmp[i]) {
 
301
    case 0: 
 
302
      _REP;
 
303
      break;
 
304
    case 1:
 
305
      _INS(1);
 
306
      break;
 
307
    case 2:
 
308
      _DEL(1);
 
309
      break;
 
310
    }
 
311
  if(X != M) _DEL(M - X)
 
312
               else if(Y != rightd) _INS(rightd-Y)
 
313
  
 
314
                                      MemFree(state[0]);
 
315
  MemFree(state);
 
316
  MemFree(tmp);
 
317
  return(v);
 
318
}
 
319
 
 
320