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