← Back to branch summary

~ubuntu-branches/ubuntu/karmic/ncbi-tools6/karmic

« back to all changes in this revision

Viewing changes to algo/blast/composition_adjustment/smith_waterman.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2006-10-22 19:32:02 UTC
  • mfrom: (1.2.1 upstream) (2.1.4 dapper)
  • Revision ID: james.westby@ubuntu.com-20061022193202-6svrvb6l52n0uhe4
Tags: 6.1.20061015-1
* New upstream release.
* Don't bother linking against indirectly used system libraries.
* Update man pages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/composition_adjustment/smith_waterman.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
 * @file smith_waterman.c
 
27
 * Routines for computing rigorous, Smith-Waterman alignments.
 
28
 *
 
29
 * @author Alejandro Schaffer, E. Michael Gertz
 
30
 */
 
31
#ifndef SKIP_DOXYGEN_PROCESSING
 
32
static char const rcsid[] =
 
33
    "$Id: smith_waterman.c,v 1.4 2006/06/29 16:50:07 gertz Exp $";
 
34
#endif /* SKIP_DOXYGEN_PROCESSING */
 
35
 
 
36
#include <algo/blast/core/ncbi_std.h>
 
37
#include <algo/blast/composition_adjustment/composition_constants.h>
 
38
#include <algo/blast/composition_adjustment/smith_waterman.h>
 
39
 
 
40
/** A structure used internally by the Smith-Waterman algorithm to
 
41
 * represent gaps */
 
42
typedef struct SwGapInfo {
 
43
    int noGap;         /**< score if opening a gap */
 
44
    int gapExists;     /**< score if continuing a gap */
 
45
} SwGapInfo;
 
46
 
 
47
 
 
48
/**
 
49
 * Compute the score and right-hand endpoints of the locally optimal
 
50
 * Smith-Waterman alignment. Called by Blast_SmithWatermanScoreOnly
 
51
 * when there are no forbidden ranges.  nonempty.  See
 
52
 * Blast_SmithWatermanScoreOnly for the meaning of the parameters to
 
53
 * this routine.
 
54
 */
 
55
static int
 
56
BLbasicSmithWatermanScoreOnly(int *score, int *matchSeqEnd, int *queryEnd,
 
57
                              const Uint1 * matchSeq, int matchSeqLength,
 
58
                              const Uint1 * query,    int queryLength,
 
59
                              int **matrix, int gapOpen, int gapExtend,
 
60
                              int positionSpecific)
 
61
{
 
62
    int bestScore;               /* best score seen so far */
 
63
    int newScore;                /* score of next entry */
 
64
    int bestMatchSeqPos, bestQueryPos; /* position ending best score in
 
65
                                          matchSeq and query sequences */
 
66
    SwGapInfo *scoreVector;      /* keeps one row of the
 
67
                                    Smith-Waterman matrix overwrite
 
68
                                    old row with new row */
 
69
    int *matrixRow;              /* one row of score matrix */
 
70
    int newGapCost;              /* cost to have a gap of one character */
 
71
    int prevScoreNoGapMatchSeq;  /* score one row and column up with
 
72
                                    no gaps */
 
73
    int prevScoreGapMatchSeq;    /* score if a gap already started in
 
74
                                    matchSeq */
 
75
    int continueGapScore;        /* score for continuing a gap in matchSeq */
 
76
    int matchSeqPos, queryPos;   /* positions in matchSeq and query */
 
77
 
 
78
    scoreVector = (SwGapInfo *) malloc(matchSeqLength * sizeof(SwGapInfo));
 
79
    if (scoreVector == NULL) {
 
80
        return -1;
 
81
    }
 
82
    bestMatchSeqPos = 0;
 
83
    bestQueryPos = 0;
 
84
    bestScore = 0;
 
85
    newGapCost = gapOpen + gapExtend;
 
86
    for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
87
        scoreVector[matchSeqPos].noGap = 0;
 
88
        scoreVector[matchSeqPos].gapExists = -gapOpen;
 
89
    }
 
90
    for (queryPos = 0;  queryPos < queryLength;  queryPos++) {
 
91
        if (positionSpecific)
 
92
            matrixRow = matrix[queryPos];
 
93
        else
 
94
            matrixRow = matrix[query[queryPos]];
 
95
        newScore = 0;
 
96
        prevScoreNoGapMatchSeq = 0;
 
97
        prevScoreGapMatchSeq = -(gapOpen);
 
98
        for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
99
            /* testing scores with a gap in matchSeq, either starting a
 
100
             * new gap or extending an existing gap*/
 
101
            if ((newScore = newScore - newGapCost) >
 
102
                (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
 
103
                prevScoreGapMatchSeq = newScore;
 
104
            /* testing scores with a gap in query, either starting a
 
105
             * new gap or extending an existing gap*/
 
106
            if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
 
107
                (continueGapScore =
 
108
                 scoreVector[matchSeqPos].gapExists - gapExtend))
 
109
                continueGapScore = newScore;
 
110
            /* compute new score extending one position in matchSeq
 
111
             * and query */
 
112
            newScore =
 
113
                prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
 
114
            if (newScore < 0)
 
115
                newScore = 0; /*Smith-Waterman locality condition*/
 
116
            /*test two alternatives*/
 
117
            if (newScore < prevScoreGapMatchSeq)
 
118
                newScore = prevScoreGapMatchSeq;
 
119
            if (newScore < continueGapScore)
 
120
                newScore = continueGapScore;
 
121
            prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap;
 
122
            scoreVector[matchSeqPos].noGap = newScore;
 
123
            scoreVector[matchSeqPos].gapExists = continueGapScore;
 
124
            if (newScore > bestScore) {
 
125
                bestScore = newScore;
 
126
                bestQueryPos = queryPos;
 
127
                bestMatchSeqPos = matchSeqPos;
 
128
            }
 
129
        }
 
130
    }
 
131
    free(scoreVector);
 
132
    if (bestScore < 0)
 
133
        bestScore = 0;
 
134
    *matchSeqEnd = bestMatchSeqPos;
 
135
    *queryEnd = bestQueryPos;
 
136
    *score = bestScore;
 
137
 
 
138
    return 0;
 
139
}
 
140
 
 
141
 
 
142
/**
 
143
 * Find the left-hand endpoints of the locally optimal Smith-Waterman
 
144
 * alignment. Called by Blast_SmithWatermanFindStart when there are no
 
145
 * forbidden ranges.  See Blast_SmithWatermanFindStartfor the meaning
 
146
 * of the parameters to this routine.
 
147
 */
 
148
static int
 
149
BLSmithWatermanFindStart(int *score_out,
 
150
                         int *matchSeqStart, int *queryStart,
 
151
                         const Uint1 * matchSeq, int matchSeqLength,
 
152
                         const Uint1 *query,
 
153
                         int **matrix, int gapOpen, int gapExtend,
 
154
                         int matchSeqEnd, int queryEnd, int score_in,
 
155
                         int positionSpecific)
 
156
{
 
157
    int bestScore;               /* best score seen so far*/
 
158
    int newScore;                /* score of next entry*/
 
159
    int bestMatchSeqPos, bestQueryPos; /*position starting best score in
 
160
                                          matchSeq and database sequences */
 
161
    SwGapInfo *scoreVector;      /* keeps one row of the Smith-Waterman
 
162
                                    matrix overwrite old row with new row */
 
163
    int *matrixRow;              /* one row of score matrix */
 
164
    int newGapCost;              /* cost to have a gap of one character */
 
165
    int prevScoreNoGapMatchSeq;  /* score one row and column up
 
166
                                    with no gaps*/
 
167
    int prevScoreGapMatchSeq;    /* score if a gap already started in
 
168
                                    matchSeq */
 
169
    int continueGapScore;        /* score for continuing a gap in query */
 
170
    int matchSeqPos, queryPos;   /* positions in matchSeq and query */
 
171
 
 
172
    scoreVector = (SwGapInfo *) malloc(matchSeqLength * sizeof(SwGapInfo));
 
173
    if (scoreVector == NULL) {
 
174
        return -1;
 
175
    }
 
176
    bestMatchSeqPos = 0;
 
177
    bestQueryPos = 0;
 
178
    bestScore = 0;
 
179
    newGapCost = gapOpen + gapExtend;
 
180
    for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
181
        scoreVector[matchSeqPos].noGap = 0;
 
182
        scoreVector[matchSeqPos].gapExists = -(gapOpen);
 
183
    }
 
184
    for (queryPos = queryEnd;  queryPos >= 0;  queryPos--) {
 
185
        if (positionSpecific)
 
186
            matrixRow = matrix[queryPos];
 
187
        else
 
188
            matrixRow = matrix[query[queryPos]];
 
189
        newScore = 0;
 
190
        prevScoreNoGapMatchSeq = 0;
 
191
        prevScoreGapMatchSeq = -(gapOpen);
 
192
        for (matchSeqPos = matchSeqEnd;  matchSeqPos >= 0;  matchSeqPos--) {
 
193
            /* testing scores with a gap in matchSeq, either starting
 
194
             * a new gap or extending an existing gap */
 
195
            if ((newScore = newScore - newGapCost) >
 
196
                (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
 
197
                prevScoreGapMatchSeq = newScore;
 
198
            /* testing scores with a gap in query, either starting a
 
199
             * new gap or extending an existing gap */
 
200
            if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
 
201
                (continueGapScore =
 
202
                 scoreVector[matchSeqPos].gapExists - gapExtend))
 
203
                continueGapScore = newScore;
 
204
            /* compute new score extending one position in matchSeq
 
205
             * and query */
 
206
            newScore =
 
207
                prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
 
208
            if (newScore < 0)
 
209
                newScore = 0; /* Smith-Waterman locality condition */
 
210
            /* test two alternatives */
 
211
            if (newScore < prevScoreGapMatchSeq)
 
212
                newScore = prevScoreGapMatchSeq;
 
213
            if (newScore < continueGapScore)
 
214
                newScore = continueGapScore;
 
215
            prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap;
 
216
            scoreVector[matchSeqPos].noGap = newScore;
 
217
            scoreVector[matchSeqPos].gapExists = continueGapScore;
 
218
            if (newScore > bestScore) {
 
219
                bestScore = newScore;
 
220
                bestQueryPos = queryPos;
 
221
                bestMatchSeqPos = matchSeqPos;
 
222
            }
 
223
            if (bestScore >= score_in)
 
224
                break;
 
225
        }
 
226
        if (bestScore >= score_in)
 
227
            break;
 
228
    }
 
229
    free(scoreVector);
 
230
    if (bestScore < 0)
 
231
        bestScore = 0;
 
232
    *matchSeqStart = bestMatchSeqPos;
 
233
    *queryStart = bestQueryPos;
 
234
    *score_out = bestScore;
 
235
 
 
236
    return 0;
 
237
}
 
238
 
 
239
 
 
240
/**
 
241
 * Compute the score and right-hand endpoints of the locally optimal
 
242
 * Smith-Waterman alignment, subject to the restriction that some
 
243
 * ranges are forbidden.  Called by Blast_SmithWatermanScoreOnly when
 
244
 * forbiddenRanges is nonempty.  See Blast_SmithWatermanScoreOnly for
 
245
 * the meaning of the parameters to this routine.
 
246
 */
 
247
static int
 
248
BLspecialSmithWatermanScoreOnly(int *score, int *matchSeqEnd, int *queryEnd,
 
249
                                const Uint1 * matchSeq, int matchSeqLength,
 
250
                                const Uint1 *query, int queryLength,
 
251
                                int **matrix, int gapOpen, int gapExtend,
 
252
                                const int *numForbidden,
 
253
                                int ** forbiddenRanges,
 
254
                                int positionSpecific)
 
255
{
 
256
    int bestScore;               /* best score seen so far */
 
257
    int newScore;                /* score of next entry*/
 
258
    int bestMatchSeqPos, bestQueryPos; /*position ending best score in
 
259
                                          matchSeq and database sequences */
 
260
    SwGapInfo *scoreVector;      /* keeps one row of the Smith-Waterman
 
261
                                    matrix overwrite old row with new row */
 
262
    int *matrixRow;              /* one row of score matrix */
 
263
    int newGapCost;              /* cost to have a gap of one character */
 
264
    int prevScoreNoGapMatchSeq;  /* score one row and column up
 
265
                                    with no gaps*/
 
266
    int prevScoreGapMatchSeq;    /* score if a gap already started in
 
267
                                    matchSeq */
 
268
    int continueGapScore;        /* score for continuing a gap in query */
 
269
    int matchSeqPos, queryPos;   /* positions in matchSeq and query */
 
270
    int forbidden;               /* is this position forbidden? */
 
271
    int f;                       /* index over forbidden positions */
 
272
 
 
273
    scoreVector = (SwGapInfo *) malloc(matchSeqLength * sizeof(SwGapInfo));
 
274
    if (scoreVector == NULL) {
 
275
        return -1;
 
276
    }
 
277
    bestMatchSeqPos = 0;
 
278
    bestQueryPos = 0;
 
279
    bestScore = 0;
 
280
    newGapCost = gapOpen + gapExtend;
 
281
    for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
282
        scoreVector[matchSeqPos].noGap = 0;
 
283
        scoreVector[matchSeqPos].gapExists = -(gapOpen);
 
284
    }
 
285
    for (queryPos = 0;  queryPos < queryLength;  queryPos++) {
 
286
        if (positionSpecific)
 
287
            matrixRow = matrix[queryPos];
 
288
        else
 
289
            matrixRow = matrix[query[queryPos]];
 
290
        newScore = 0;
 
291
        prevScoreNoGapMatchSeq = 0;
 
292
        prevScoreGapMatchSeq = -(gapOpen);
 
293
        for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
294
            /* testing scores with a gap in matchSeq, either starting
 
295
             * a new gap or extending an existing gap */
 
296
            if ((newScore = newScore - newGapCost) >
 
297
                (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
 
298
                prevScoreGapMatchSeq = newScore;
 
299
            /* testing scores with a gap in query, either starting a
 
300
             * new gap or extending an existing gap */
 
301
            if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
 
302
                (continueGapScore =
 
303
                 scoreVector[matchSeqPos].gapExists - gapExtend))
 
304
                continueGapScore = newScore;
 
305
            /* compute new score extending one position in matchSeq
 
306
             * and query */
 
307
            forbidden = FALSE;
 
308
            for (f = 0;  f < numForbidden[queryPos];  f++) {
 
309
                if ((matchSeqPos >= forbiddenRanges[queryPos][2 * f]) &&
 
310
                    (matchSeqPos <= forbiddenRanges[queryPos][2*f + 1])) {
 
311
                    forbidden = TRUE;
 
312
                    break;
 
313
                }
 
314
            }
 
315
            if (forbidden)
 
316
                newScore = COMPO_SCORE_MIN;
 
317
            else
 
318
                newScore =
 
319
                    prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
 
320
            if (newScore < 0)
 
321
                newScore = 0; /* Smith-Waterman locality condition */
 
322
            /* test two alternatives */
 
323
            if (newScore < prevScoreGapMatchSeq)
 
324
                newScore = prevScoreGapMatchSeq;
 
325
            if (newScore < continueGapScore)
 
326
                newScore = continueGapScore;
 
327
            prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap;
 
328
            scoreVector[matchSeqPos].noGap = newScore;
 
329
            scoreVector[matchSeqPos].gapExists = continueGapScore;
 
330
            if (newScore > bestScore) {
 
331
                bestScore = newScore;
 
332
                bestQueryPos = queryPos;
 
333
                bestMatchSeqPos = matchSeqPos;
 
334
            }
 
335
        }
 
336
    }
 
337
    free(scoreVector);
 
338
    if (bestScore < 0)
 
339
        bestScore = 0;
 
340
    *matchSeqEnd = bestMatchSeqPos;
 
341
    *queryEnd = bestQueryPos;
 
342
    *score = bestScore;
 
343
 
 
344
    return 0;
 
345
}
 
346
 
 
347
 
 
348
/**
 
349
 * Find the left-hand endpoints of the locally optimal Smith-Waterman
 
350
 * alignment, subject to the restriction that certain ranges may not
 
351
 * be aligned. Called by Blast_SmithWatermanFindStart if
 
352
 * forbiddenRanges is nonempty.  See Blast_SmithWatermanFindStartfor
 
353
 * the meaning of the parameters to this routine.
 
354
 */
 
355
static int
 
356
BLspecialSmithWatermanFindStart(int * score_out,
 
357
                                int *matchSeqStart, int *queryStart,
 
358
                                const Uint1 * matchSeq, int matchSeqLength,
 
359
                                const Uint1 *query, int **matrix,
 
360
                                int gapOpen, int gapExtend, int matchSeqEnd,
 
361
                                int queryEnd, int score_in,
 
362
                                const int *numForbidden,
 
363
                                int ** forbiddenRanges,
 
364
                                int positionSpecific)
 
365
{
 
366
    int bestScore;               /* best score seen so far */
 
367
    int newScore;                /* score of next entry */
 
368
    int bestMatchSeqPos, bestQueryPos; /* position starting best score in
 
369
                                          matchSeq and database sequences */
 
370
    SwGapInfo *scoreVector;      /* keeps one row of the
 
371
                                    Smith-Waterman matrix; overwrite
 
372
                                    old row with new row*/
 
373
    int *matrixRow;              /* one row of score matrix */
 
374
    int newGapCost;              /* cost to have a gap of one character */
 
375
    int prevScoreNoGapMatchSeq;  /* score one row and column up
 
376
                                    with no gaps*/
 
377
    int prevScoreGapMatchSeq;    /* score if a gap already started in
 
378
                                    matchSeq */
 
379
    int continueGapScore;        /* score for continuing a gap in query */
 
380
    int matchSeqPos, queryPos;   /* positions in matchSeq and query */
 
381
    int forbidden;               /* is this position forbidden? */
 
382
    int f;                       /* index over forbidden positions */
 
383
    
 
384
    scoreVector = (SwGapInfo *) malloc(matchSeqLength * sizeof(SwGapInfo));
 
385
    if (scoreVector == NULL) {
 
386
        return -1;
 
387
    }
 
388
    bestMatchSeqPos = 0;
 
389
    bestQueryPos = 0;
 
390
    bestScore = 0;
 
391
    newGapCost = gapOpen + gapExtend;
 
392
    for (matchSeqPos = 0;  matchSeqPos < matchSeqLength;  matchSeqPos++) {
 
393
        scoreVector[matchSeqPos].noGap = 0;
 
394
        scoreVector[matchSeqPos].gapExists = -(gapOpen);
 
395
    }
 
396
    for (queryPos = queryEnd;  queryPos >= 0;  queryPos--) {
 
397
        if (positionSpecific)
 
398
            matrixRow = matrix[queryPos];
 
399
        else
 
400
            matrixRow = matrix[query[queryPos]];
 
401
        newScore = 0;
 
402
        prevScoreNoGapMatchSeq = 0;
 
403
        prevScoreGapMatchSeq = -(gapOpen);
 
404
        for (matchSeqPos = matchSeqEnd;  matchSeqPos >= 0;  matchSeqPos--) {
 
405
            /* testing scores with a gap in matchSeq, either starting a
 
406
             * new gap or extending an existing gap*/
 
407
            if ((newScore = newScore - newGapCost) >
 
408
                (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
 
409
                prevScoreGapMatchSeq = newScore;
 
410
            /* testing scores with a gap in query, either starting a
 
411
             * new gap or extending an existing gap*/
 
412
            if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
 
413
                (continueGapScore =
 
414
                 scoreVector[matchSeqPos].gapExists - gapExtend))
 
415
                continueGapScore = newScore;
 
416
            /* compute new score extending one position in matchSeq
 
417
             * and query */
 
418
            forbidden = FALSE;
 
419
            for (f = 0;  f < numForbidden[queryPos];  f++) {
 
420
                if ((matchSeqPos >= forbiddenRanges[queryPos][2 * f]) &&
 
421
                    (matchSeqPos <= forbiddenRanges[queryPos][2*f + 1])) {
 
422
                    forbidden = TRUE;
 
423
                    break;
 
424
                }
 
425
            }
 
426
            if (forbidden)
 
427
                newScore = COMPO_SCORE_MIN;
 
428
            else
 
429
                newScore =
 
430
                    prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
 
431
            if (newScore < 0)
 
432
                newScore = 0; /* Smith-Waterman locality condition */
 
433
            /* test two alternatives */
 
434
            if (newScore < prevScoreGapMatchSeq)
 
435
                newScore = prevScoreGapMatchSeq;
 
436
            if (newScore < continueGapScore)
 
437
                newScore = continueGapScore;
 
438
            prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap;
 
439
            scoreVector[matchSeqPos].noGap = newScore;
 
440
            scoreVector[matchSeqPos].gapExists = continueGapScore;
 
441
            if (newScore > bestScore) {
 
442
                bestScore = newScore;
 
443
                bestQueryPos = queryPos;
 
444
                bestMatchSeqPos = matchSeqPos;
 
445
            }
 
446
            if (bestScore >= score_in)
 
447
                break;
 
448
        }
 
449
        if (bestScore >= score_in)
 
450
            break;
 
451
    }
 
452
    free(scoreVector);
 
453
    if (bestScore < 0)
 
454
        bestScore = 0;
 
455
    *matchSeqStart = bestMatchSeqPos;
 
456
    *queryStart = bestQueryPos;
 
457
    *score_out = bestScore;
 
458
    
 
459
    return 0;
 
460
}
 
461
 
 
462
 
 
463
/* Documented in smith_waterman.h. */
 
464
void
 
465
Blast_ForbiddenRangesRelease(Blast_ForbiddenRanges * self)
 
466
{
 
467
    int f;
 
468
    if (self->ranges) {
 
469
        for (f = 0;  f < self->capacity;  f++) free(self->ranges[f]);
 
470
    }
 
471
    free(self->ranges);       self->ranges       = NULL;
 
472
    free(self->numForbidden); self->numForbidden = NULL;
 
473
}
 
474
 
 
475
 
 
476
/* Documented in smith_waterman.h. */
 
477
int
 
478
Blast_ForbiddenRangesInitialize(Blast_ForbiddenRanges * self,
 
479
                                int capacity)
 
480
{
 
481
    int f;
 
482
    self->capacity  = capacity;
 
483
    self->numForbidden = NULL;
 
484
    self->ranges       = NULL;
 
485
    self->isEmpty      = TRUE;
 
486
 
 
487
    self->numForbidden = (int *) calloc(capacity, sizeof(int));
 
488
    if (self->numForbidden == NULL)
 
489
        goto error_return;
 
490
    self->ranges       = (int **) calloc(capacity, sizeof(int *));
 
491
    if (self->ranges == NULL) 
 
492
        goto error_return;
 
493
    for (f = 0;  f < capacity;  f++) {
 
494
        self->numForbidden[f] = 0;
 
495
        self->ranges[f]       = (int *) malloc(2 * sizeof(int));
 
496
        if (self->ranges[f] == NULL) 
 
497
            goto error_return;
 
498
        self->ranges[f][0]    = 0;
 
499
        self->ranges[f][1]    = 0;
 
500
    }
 
501
    return 0;
 
502
error_return:
 
503
    Blast_ForbiddenRangesRelease(self);
 
504
    return -1;
 
505
}
 
506
 
 
507
 
 
508
/* Documented in smith_waterman.h. */
 
509
void Blast_ForbiddenRangesClear(Blast_ForbiddenRanges * self)
 
510
{
 
511
    int f;
 
512
    for (f = 0;  f < self->capacity;  f++) {
 
513
        self->numForbidden[f] = 0;
 
514
    }
 
515
    self->isEmpty = TRUE;
 
516
}
 
517
 
 
518
 
 
519
/* Documented in smith_waterman.h. */
 
520
int
 
521
Blast_ForbiddenRangesPush(Blast_ForbiddenRanges * self,
 
522
                          int queryStart,
 
523
                          int queryEnd,
 
524
                          int matchStart,
 
525
                          int matchEnd)
 
526
{
 
527
    int f;
 
528
    for (f = queryStart;  f < queryEnd;  f++) {
 
529
        int last = 2 * self->numForbidden[f];
 
530
        if (0 != last) {    /* we must resize the array */
 
531
            int * new_ranges =
 
532
                realloc(self->ranges[f], (last + 2) * sizeof(int));
 
533
            if (new_ranges == NULL) 
 
534
                return -1;
 
535
            self->ranges[f] = new_ranges;
 
536
        }
 
537
        self->ranges[f][last]     = matchStart;
 
538
        self->ranges[f][last + 1] = matchEnd;
 
539
 
 
540
        self->numForbidden[f]++;
 
541
    }
 
542
    self->isEmpty = FALSE;
 
543
 
 
544
    return 0;
 
545
}
 
546
 
 
547
 
 
548
/* Documented in smith_waterman.h. */
 
549
int
 
550
Blast_SmithWatermanScoreOnly(int *score,
 
551
                             int *matchSeqEnd, int *queryEnd,
 
552
                             const Uint1 * subject_data, int subject_length,
 
553
                             const Uint1 * query_data, int query_length,
 
554
                             int **matrix,
 
555
                             int gapOpen,
 
556
                             int gapExtend,
 
557
                             int positionSpecific,
 
558
                             const Blast_ForbiddenRanges * forbiddenRanges )
 
559
{
 
560
    if (forbiddenRanges->isEmpty) {
 
561
        return BLbasicSmithWatermanScoreOnly(score, matchSeqEnd,
 
562
                                             queryEnd, subject_data,
 
563
                                             subject_length,
 
564
                                             query_data, query_length,
 
565
                                             matrix, gapOpen,
 
566
                                             gapExtend,
 
567
                                             positionSpecific);
 
568
    } else {
 
569
        return BLspecialSmithWatermanScoreOnly(score, matchSeqEnd,
 
570
                                               queryEnd, subject_data,
 
571
                                               subject_length,
 
572
                                               query_data,
 
573
                                               query_length, matrix,
 
574
                                               gapOpen, gapExtend,
 
575
                                               forbiddenRanges->numForbidden,
 
576
                                               forbiddenRanges->ranges,
 
577
                                               positionSpecific);
 
578
    }
 
579
}
 
580
 
 
581
 
 
582
/* Documented in smith_waterman.h. */
 
583
int
 
584
Blast_SmithWatermanFindStart(int * score_out,
 
585
                             int *matchSeqStart,
 
586
                             int *queryStart,
 
587
                             const Uint1 * subject_data, int subject_length,
 
588
                             const Uint1 * query_data,
 
589
                             int **matrix,
 
590
                             int gapOpen,
 
591
                             int gapExtend,
 
592
                             int matchSeqEnd,
 
593
                             int queryEnd,
 
594
                             int score_in,
 
595
                             int positionSpecific,
 
596
                             const Blast_ForbiddenRanges * forbiddenRanges)
 
597
{
 
598
    if (forbiddenRanges->isEmpty) {
 
599
        return BLSmithWatermanFindStart(score_out, matchSeqStart,
 
600
                                        queryStart, subject_data,
 
601
                                        subject_length, query_data,
 
602
                                        matrix, gapOpen, gapExtend,
 
603
                                        matchSeqEnd, queryEnd,
 
604
                                        score_in, positionSpecific);
 
605
    } else {
 
606
        return BLspecialSmithWatermanFindStart(score_out,
 
607
                                               matchSeqStart,
 
608
                                               queryStart,
 
609
                                               subject_data,
 
610
                                               subject_length,
 
611
                                               query_data, matrix,
 
612
                                               gapOpen, gapExtend,
 
613
                                               matchSeqEnd, queryEnd,
 
614
                                               score_in,
 
615
                                               forbiddenRanges->numForbidden,
 
616
                                               forbiddenRanges->ranges,
 
617
                                               positionSpecific);
 
618
    }
 
619
}