1 static char const rcsid[] = "$Id: bandalgn.c,v 6.17 2003/05/30 17:25:36 coulouri Exp $";
2
3 /* ===========================================================================
4 *
5 * PUBLIC DOMAIN NOTICE
6 * National Center for Biotechnology Information
7 *
8 * This software/database is a "United States Government Work" under the
9 * terms of the United States Copyright Act. It was written as part of
10 * the author's official duties as a United States Government employee and
11 * thus cannot be copyrighted. This software/database is freely available
12 * to the public for use. The National Library of Medicine and the U.S.
13 * Government have not placed any restriction on its use or reproduction.
14 *
15 * Although all reasonable efforts have been taken to ensure the accuracy
16 * and reliability of the software and data, the NLM and the U.S.
17 * Government do not and cannot warrant the performance or results that
18 * may be obtained by using this software or data. The NLM and the U.S.
19 * Government disclaim all warranties, express or implied, including
20 * warranties of performance, merchantability or fitness for any particular
21 * purpose.
22 *
23 * Please cite the author in any work or product based on this material.
24 *
25 * ===========================================================================*/
26 /*****************************************************************************
27
28 File name: bandalgn.c
29
30 Author: Gennadiy Savchuk, Jinqhui Zhang, Tom Madden
31
32 Contents: Functions to perform both local and global banded alignments.
33
34 ****************************************************************************/
35 /* $Revision: 6.17 $
36 *
37 * $Log: bandalgn.c,v $
38 * Revision 6.17 2003/05/30 17:25:36 coulouri
39 * add rcsid
40 *
41 * Revision 6.16 2001/05/25 19:32:48 vakatov
42 * Nested comment typo fixed
43 *
44 * Revision 6.15 2000/02/10 22:47:06 vakatov
45 * DLL'zation for MSVC on PC, Win-NT
46 *
47 * Revision 6.14 1999/09/01 21:46:49 chappey
48 * SeqPortFree in load_data
49 *
50 * Revision 6.13 1999/07/13 15:39:48 sicotte
51 * Fixed casting warnings
52 *
53 * Revision 6.12 1999/04/06 15:21:00 sicotte
54 * Fixed prototype problems and lint casting and.
55 *
56 * Revision 6.11 1998/08/24 20:20:32 kans
57 * fixed -v -fd warnings
58 *
59 * Revision 6.10 1998/05/04 19:58:21 tatiana
60 * change SEV_ERROR to SEV_WARNING
61 *
62 * Revision 6.9 1998/01/16 21:29:43 chappey
63 * Remove function CC_GetExtremes and use now SeqAlignStart, SeqAlignStop in salsap.c
64 *
65 * Revision 6.8 1998/01/02 21:58:55 chappey
66 * printf removed
67 *
68 * Revision 6.7 1997/12/19 20:47:16 chappey
69 * bugs fixed
70 *
71 * Revision 6.6 1997/10/21 23:09:39 chappey
72 * fixes on ChangeGlobalBandMatrix by Colombe and Hugues
73 *
74 * Revision 6.4 1997/10/14 16:21:40 tatiana
75 * merge utility functions
76 *
77 * Revision 6.3 1997/10/02 15:54:54 kans
78 * added two includes and a return value
79 *
80 * Revision 6.2 1997/10/02 15:17:55 tatiana
81 * global align utility added
82 *
83 * Revision 6.1 1997/09/16 19:29:14 kans
84 * minor fix (CC)
85 *
86 * Revision 6.0 1997/08/25 18:52:12 madden
87 * Revision changed to 6.0
88 *
89 * Revision 1.3 1997/06/23 16:15:52 tatiana
90 * GlobalBandAlign struct changed to use SeqLocs instead of SeqIds
91 *
92 * Revision 1.2 1997/03/05 17:31:21 savchuk
93 * Sliced for 5 pieces (bandalg[n,0-5].c)
94 * to fit the Windoze MSVC 16-bit compiler limitations.
95 *
96 * Revision 1.1 1997/01/22 14:11:05 madden
97 * Initial revision
98 *
99 */
100
101
102
103 /* A PACKAGE FOR GLOBALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
104
105 To invoke, call ALIGN(A,B,M,N,L,U,W,G,H,S).
106 The parameters are explained as follows:
107 A, B : two sequences to be aligned
108 M : the length of sequence A
109 N : the length of sequence B
110 L : lower bound of the band
111 U : upper bound of the band
112 W : scoring table for matches and mismatches
113 G : gap-opening penalty
114 H : gap-extension penalty
115 S : script for DISPLAY routine
116 */
117
118 #include <seqport.h>
119 #include <blast.h>
120 #include <salutil.h>
121 #include <bandalgn.h>
122
123
124 NLM_EXTERN PSUGapOptionsPtr
PSUGapOptionsDelete(PSUGapOptionsPtr options)125 PSUGapOptionsDelete(PSUGapOptionsPtr options)
126
127 {
128 options = MemFree(options);
129
130 return NULL;
131 }
132
133 /*
134 Create the PSUGapOptionsPtr and fill in with default values.
135 */
136
137 NLM_EXTERN PSUGapOptionsPtr
PSUGapOptionsCreate(Uint1 search_type)138 PSUGapOptionsCreate(Uint1 search_type)
139
140 {
141 PSUGapOptionsPtr options;
142
143 options = (PSUGapOptionsPtr) MemNew(sizeof(PSUGapOptions));
144 if (options == NULL)
145 return NULL;
146
147 /* defaults from WEbb lab. */
148 /* Should be defines above?? */
149 options->up = 282;
150 options->low = -2329;
151 options->gopen = 5;
152 options->gext = 1;
153 options->gext = 1;
154 options->start_diag = options->low;
155 options->width = options->up - options->low + 1;
156
157 return options;
158 }
159
160 /*
161 Deletes the GlobalBandStruct, including the options.
162 Does not delete the sequence or the ID's.
163 */
164
165
166 NLM_EXTERN GlobalBandStructPtr
GlobalBandStructDelete(GlobalBandStructPtr gbsp)167 GlobalBandStructDelete(GlobalBandStructPtr gbsp)
168
169 {
170 if (gbsp == NULL)
171 return NULL;
172
173 gbsp->options = PSUGapOptionsDelete(gbsp->options);
174
175 gbsp->edit_block = GapXEditBlockDelete(gbsp->edit_block);
176
177 gbsp = MemFree(gbsp);
178
179 return NULL;
180 }
181
182 /*
183 Create the GlobalBandStructPtr and defaults for a search_type.
184 */
185
186 NLM_EXTERN GlobalBandStructPtr
GlobalBandStructCreate(Uint1 search_type)187 GlobalBandStructCreate(Uint1 search_type)
188
189 {
190 GlobalBandStructPtr gbsp;
191
192 gbsp = (GlobalBandStructPtr) MemNew(sizeof(GlobalBandStruct));
193 if (gbsp == NULL)
194 return NULL;
195
196 gbsp->search_type = search_type;
197 gbsp->options = PSUGapOptionsCreate(search_type);
198 if (gbsp->options == NULL)
199 {
200 gbsp = GlobalBandStructDelete(gbsp);
201 return NULL;
202 }
203
204 return gbsp;
205 }
206
207 /*
208 Performs a global alignment, producing a SeqAlign.
209 */
210 NLM_EXTERN SeqAlignPtr
GlobalBandToSeqAlign(GlobalBandStructPtr gbsp)211 GlobalBandToSeqAlign (GlobalBandStructPtr gbsp)
212 {
213
214 SeqAlignPtr seqalign;
215 GapXEditBlockPtr edit_block;
216 SeqLocPtr l1, l2;
217
218 if (gbsp == NULL)
219 return NULL;
220
221 if (GlobalBandToEditBlock(gbsp) == FALSE)
222 {
223 return NULL;
224 }
225
226 if (gbsp->edit_block == NULL)
227 return NULL;
228 edit_block = gbsp->edit_block;
229 l1 = gbsp->seqloc1;
230 l2 = gbsp->seqloc2;
231 edit_block->length1 = SeqLocLen(l1);
232 edit_block->length2 = SeqLocLen(l2);
233 if (SeqLocStrand(l1) == Seq_strand_minus) {
234 edit_block->frame1 = -1;
235 } else {
236 edit_block->frame1 = 1;
237 }
238 if (SeqLocStrand(l2) == Seq_strand_minus) {
239 edit_block->frame2 = -1;
240 } else {
241 edit_block->frame2 = 1;
242 }
243 seqalign = GapXEditBlockToSeqAlign(gbsp->edit_block, SeqLocId(l2), SeqLocId(l1));
244
245 gbsp->edit_block = GapXEditBlockDelete(gbsp->edit_block);
246
247 return seqalign;
248 }
249
250 /*
251 Performs a global alignment, producing an EditBlock, which
252 can be made into a SeqAlign.
253 */
254
255 NLM_EXTERN Boolean
GlobalBandToEditBlock(GlobalBandStructPtr gbsp)256 GlobalBandToEditBlock(GlobalBandStructPtr gbsp)
257 {
258 GapXEditBlockPtr edit_block;
259 FloatHi score;
260 Int4 Slen;
261 Int4Ptr S; /* conversion operations */
262
263 if (gbsp == NULL)
264 return FALSE;
265
266 S = MemGet((gbsp->seq1_length+gbsp->seq2_length)*sizeof(Int4), MGET_ERRPOST);
267
268 switch(gbsp->search_type)
269 {
270
271 case G_BAND_LINEAR:
272 score = (FloatHi) gband_linear(gbsp->seq1-1, gbsp->seq2-1,
273 gbsp->seq1_length, gbsp->seq2_length,
274 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
275 break;
276 case G_BAND_QUADRATIC:
277 score = (FloatHi) gband_quadratic(gbsp->seq1-1, gbsp->seq2-1,
278 gbsp->seq1_length, gbsp->seq2_length,
279 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
280 break;
281 case G_BAND_LGAP:
282 score = (FloatHi) gband_linear_gap(gbsp->seq1-1, gbsp->seq2-1,
283 gbsp->seq1_length, gbsp->seq2_length,
284 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
285 break;
286 case G_BAND_QGAP:
287 score = (FloatHi) gband_linear_qgap(gbsp->seq1-1, gbsp->seq2-1,
288 gbsp->seq1_length, gbsp->seq2_length,
289 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
290 break;
291 case G_BAND_L3GAP:
292 score = (FloatHi) gband_l3gap(gbsp->seq1-1, gbsp->seq2-1,
293 gbsp->seq1_length, gbsp->seq2_length,
294 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
295 break;
296 case G_BAND_Q3GAP:
297 score = (FloatHi) gband_q3gap(gbsp->seq1-1, gbsp->seq2-1,
298 gbsp->seq1_length, gbsp->seq2_length,
299 gbsp->matrix, gbsp->options, S, &Slen)/BND_DIGIT;
300 break;
301 default:
302 score = 0;
303 ErrPostEx(SEV_ERROR, 0, 0, "Unknown method.");
304 return FALSE;
305 }
306 if (score == 0) {
307 return FALSE;
308 }
309 edit_block = TracebackToGapXEditBlock(gbsp->seq1-1, gbsp->seq2-1, gbsp->seq1_length, gbsp->seq2_length, S, 0, 0);
310
311 S = MemFree(S);
312
313 gbsp->edit_block = edit_block;
314 gbsp->score = (Int4) score;
315 gbsp->alignment_length = Slen;
316
317 return TRUE;
318 }
319
320 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
321
322 To invoke, call BAND_LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,&SI,&SJ,&EI,&EJ).
323 The parameters are explained as follows:
324 A, B : two sequences to be aligned
325 M : the length of sequence A
326 N : the length of sequence B
327 L : lower bound of the band
328 U : upper bound of the band
329 W : scoring table for matches and mismatches
330 G : gap-opening penalty
331 H : gap-extension penalty
332 S : script for DISPLAY routine
333 *SI : starting position of sequence A in the optimal local alignment
334 *SJ : starting position of sequence B in the optimal local alignment
335 *EI : ending position of sequence A in the optimal local alignment
336 *EJ : ending position of sequence B in the optimal local alignment
337 */
338
339 /*
340 Deletes the LocalBandStruct, including the options.
341 Does not delete the sequence or the ID's.
342 */
343
344 NLM_EXTERN LocalBandStructPtr
LocalBandStructDelete(LocalBandStructPtr lbsp)345 LocalBandStructDelete(LocalBandStructPtr lbsp)
346
347 {
348 if (lbsp == NULL)
349 return NULL;
350
351 lbsp->options = PSUGapOptionsDelete(lbsp->options);
352
353 lbsp->edit_block = GapXEditBlockDelete(lbsp->edit_block);
354
355 lbsp = MemFree(lbsp);
356
357 return NULL;
358 }
359
360 /*
361 Create the LocalBandStructPtr and defaults for a search_type.
362 */
363
364 NLM_EXTERN LocalBandStructPtr
LocalBandStructCreate(Uint1 search_type)365 LocalBandStructCreate(Uint1 search_type)
366
367 {
368 LocalBandStructPtr lbsp;
369
370 lbsp = (LocalBandStructPtr) MemNew(sizeof(LocalBandStruct));
371 if (lbsp == NULL)
372 return NULL;
373
374 lbsp->search_type = search_type;
375 lbsp->options = PSUGapOptionsCreate(search_type);
376 if (lbsp->options == NULL)
377 {
378 lbsp = LocalBandStructDelete(lbsp);
379 return NULL;
380 }
381
382 return lbsp;
383 }
384
385 /*
386 Performs local alignment, producing a SeqAlign.
387 */
388 NLM_EXTERN SeqAlignPtr
LocalBandToSeqAlign(LocalBandStructPtr lbsp)389 LocalBandToSeqAlign (LocalBandStructPtr lbsp)
390 {
391
392 SeqAlignPtr seqalign;
393
394 GapXEditBlockPtr edit_block;
395 SeqLocPtr l1, l2;
396
397 if (lbsp == NULL)
398 return NULL;
399
400 if (LocalBandToEditBlock(lbsp) == FALSE)
401 {
402 return NULL;
403 }
404
405 if (lbsp->edit_block == NULL)
406 return NULL;
407 edit_block = lbsp->edit_block;
408 l1 = lbsp->seqloc1;
409 l2 = lbsp->seqloc2;
410 if (SeqLocStrand(l1) == Seq_strand_minus) {
411 edit_block->frame1 = -1;
412 } else {
413 edit_block->frame1 = 1;
414 }
415 if (SeqLocStrand(l2) == Seq_strand_minus) {
416 edit_block->frame2 = -1;
417 } else {
418 edit_block->frame2 = 1;
419 }
420 seqalign = GapXEditBlockToSeqAlign(lbsp->edit_block, SeqLocId(l2), SeqLocId(l1));
421
422 lbsp->edit_block = GapXEditBlockDelete(lbsp->edit_block);
423
424 return seqalign;
425 }
426
427 /*
428 Performs a global alignment, producing an EditBlock, which
429 can be made into a SeqAlign.
430 */
431
432 NLM_EXTERN Boolean
LocalBandToEditBlock(LocalBandStructPtr lbsp)433 LocalBandToEditBlock(LocalBandStructPtr lbsp)
434 {
435 GapXEditBlockPtr edit_block;
436 FloatHi score;
437 Int4Ptr S; /* conversion operations */
438 Uint1Ptr seq1, seq2;
439
440 if (lbsp == NULL)
441 return FALSE;
442
443 S = (Int4Ptr)MemGet((lbsp->seq1_length + lbsp->seq2_length) * sizeof(Int4), MGET_ERRPOST);
444 score = BAND_LOCAL_ALIGN(lbsp->seq1-1, lbsp->seq2-1, lbsp->seq1_length, lbsp->seq2_length, lbsp->matrix, lbsp->options, S, &(lbsp->seq1_start), &(lbsp->seq2_start), &(lbsp->seq1_end), &(lbsp->seq2_end), lbsp->search_type);
445
446 /* Compensate for one-offset. */
447 lbsp->seq1_start--;
448 lbsp->seq2_start--;
449 lbsp->seq1_end--;
450 lbsp->seq2_end--;
451
452 seq1 = lbsp->seq1;
453 seq2 = lbsp->seq2;
454
455 seq1 += lbsp->seq1_start;
456 seq2 += lbsp->seq2_start;
457
458 if (score > 0)
459 edit_block = TracebackToGapXEditBlock(seq1-1, seq2-1, (lbsp->seq1_length-lbsp->seq1_start), (lbsp->seq2_length-lbsp->seq2_start), S, lbsp->seq1_start, lbsp->seq2_start);
460
461 S = MemFree(S);
462
463 lbsp->edit_block = edit_block;
464 lbsp->score = (Int4)score;
465
466 return TRUE;
467 }
468
469 Int4
BAND_LOCAL_ALIGN(Uint1Ptr A,Uint1Ptr B,Int4 M,Int4 N,Int4Ptr PNTR matrix,PSUGapOptionsPtr options,Int4Ptr S,Int4Ptr psi,Int4Ptr psj,Int4Ptr pei,Int4Ptr pej,Int4 align_type)470 BAND_LOCAL_ALIGN(Uint1Ptr A, Uint1Ptr B,
471 Int4 M, Int4 N,
472 Int4Ptr PNTR matrix, PSUGapOptionsPtr options,
473 Int4Ptr S,
474 Int4Ptr psi, Int4Ptr psj,
475 Int4Ptr pei, Int4Ptr pej,
476 Int4 align_type)
477 {
478 Int4 band;
479 Boolean flag;
480 register Int4 i, j, si, ei, ib;
481 register Int4 c, d, e, t, m;
482 Int4 leftd, rightd, G, H, low, up;
483 Int4 best_score, score1;
484 Int4 starti = 1, startj = 1, endi = M - 1, endj = N - 1;
485 register Int4Ptr wa;
486 register Int4 curd;
487 register dp_ptr dp;
488 Int4 Slen;
489 data_t data;
490 Int4Ptr PNTR W;
491
492 G = options->gopen;
493 H = options->gext;
494 m = G + H;
495 low = MAX(-M, options->low);
496 up = MIN(N, options->up);
497 W = matrix;
498
499 if (N <= 0) {
500 *psi = *psj = *pei = *pej;
501 return -1;
502 }
503 if (M <= 0) {
504 *psi = *psj = *pei = *pej;
505 return -1;
506 }
507 band = up - low + 1;
508 if (band < 1) {
509 ErrPostEx(SEV_WARNING, 0, 0, "low > up is unacceptable");
510 return -1;
511 }
512 j = (band+2) * sizeof(dp_node);
513 data.CD = (dp_ptr) MemGet(j, MGET_ERRPOST);
514
515 if (low > 0) leftd = 1;
516 else if (up < 0) leftd = band;
517 else leftd = 1-low;
518 rightd = band;
519 si = MAX(0, -up);
520 ei = MIN(M, N - low);
521 data.CD[leftd].CC = 0;
522 for (j = leftd+1; j <= rightd; j++) {
523 data.CD[j].CC = 0;
524 data.CD[j].DD = -m;
525 }
526 data.CD[rightd+1].CC = MININT;
527 data.CD[rightd].DD = MININT;
528 best_score = 0;
529 endi = si;
530 endj = si + low;
531 data.CD[leftd-1].CC = MININT;
532 data.CD[leftd].DD = -m; c = 0;
533 for (i = si+1; i <= ei; i++) {
534 score1 = best_score;
535 if (i > N-up) rightd--;
536 if (leftd > 1) leftd--;
537 wa = W[A[i]];
538 d = data.CD[leftd].DD;
539 if ((ib = leftd + low - 1 + i) > 0) c = data.CD[leftd].CC+wa[B[ib]];
540 if (c < 0) c = 0;
541 else if (c > best_score) {
542 best_score = c;
543 endi = i;
544 endj = ib;
545 }
546 e = c-m;
547 data.CD[leftd].CC = c;
548 for (curd=leftd+1, dp = &data.CD[curd], curd += i+low-1;
549 curd <= rightd+i+low-1; curd++) {
550 c = dp->CC + wa[B[curd]];
551 if ((d=dp->DD) > c) c = d;
552 if (e > c) {
553 if (e > 0) {
554 dp->CC = e; e-=H;
555 *(((Int4Ptr) (dp++))-1) = d-H;
556 } else {
557 dp->CC = *(((Int4Ptr) (dp))-1) = 0;
558 dp++;
559 }
560 } else {
561 if (c <= 0) {
562 dp->CC = *(((Int4Ptr) (dp))-1) = 0;
563 dp++;
564 } else {
565 if (c > best_score) {
566 best_score = c;
567 endj = curd;
568 }
569 dp->CC = c;
570 if ((c -= m) > (e -= H)) e= c;
571 if (c > (d-=H))
572 *((Int4Ptr) (dp++)-1) = c;
573 else *((Int4Ptr) (dp++)-1) = d;
574 }
575 }
576 }
577 if (score1 < best_score) endi = i;
578 }
579 leftd = MAX(1, -endi - low + 1);
580 rightd = band - (up - (endj - endi));
581 data.CD[rightd].CC = 0; data.CD[rightd + 1].DD = -m;
582 t = -G;
583 for (j = rightd - 1; j >= leftd; j--) {
584 data.CD[j].CC = t = t - H;
585 data.CD[j + 1].DD = t - m;
586 }
587 for (j = rightd + 1; j <= band; ++j) data.CD[j].CC = MININT;
588 data.CD[leftd - 1].CC = data.CD[leftd].DD = MININT;
589 flag = FALSE;
590 for (i = endi; i >= 1; i--) {
591 if (i + low <= 0) leftd++;
592 if (rightd < band) rightd++;
593 wa = W[A[i]];
594 d = data.CD[rightd].DD;
595 if ((ib = rightd + low - 1 + i) <= N) c = data.CD[rightd].CC + wa[B[ib]];
596 else c = MININT;
597 if (d > c) c = d;
598 e = c - m;
599 data.CD[rightd].CC = c;
600 if ((d -= H) < e) data.CD[rightd + 1].DD = e;
601 else data.CD[rightd + 1].DD = d;
602 if (c == best_score) {
603 starti = i;
604 startj = ib;
605 flag = TRUE;
606 break;
607 }
608 for (curd = rightd - 1, dp = &data.CD[curd]; curd >= leftd; curd--) {
609 c = dp->CC + wa[B[curd + low + i - 1]];
610 if ((d = dp->DD) > c) c = d;
611 if (e >= c) {
612 dp->CC = e; e -= H;
613 (dp-- + 1)->DD = d - H;
614 continue;
615 }
616 dp->CC = c;
617 if (c == best_score) {
618 starti = i;
619 startj = curd + low + i - 1;
620 flag = TRUE;
621 break;
622 }
623 if ((c -= m) > (e -= H)) e = c;
624 if (c < (d -= H)) (dp-- + 1)->DD = d; else (dp-- + 1)->DD = c;
625 }
626 if (flag == TRUE)
627 break;
628 }
629 MemFree(data.CD);
630 if (starti < 0 || starti > M || startj < 0 || startj > N)
631 {
632 ErrPostEx(SEV_WARNING, 0, 0, "starti=%d, startj=%d\n", starti, startj);
633 *psi = *psj = *pei = *pej;
634 return -1;
635 }
636 *psi = starti;
637 *psj = startj;
638 *pei = endi;
639 *pej = endj;
640
641 options->start_diag = low - (startj - starti);
642 options->width = up - (startj - starti) - low + 1;
643
644 switch(align_type) {
645 /* XXX Float Score is converted to Int4 --> Round-off Errors */
646 case L_BAND_LINEAR:
647 return (Int4) (gband_linear(A + starti - 1, B + startj - 1,
648 endi - starti + 1, endj - startj + 1,
649 matrix, options, S, &Slen) / BND_DIGIT);
650 case L_BAND_QUADRATIC:
651 return (Int4) ( gband_quadratic(A + starti - 1, B + startj - 1,
652 endi - starti + 1, endj - startj + 1,
653 matrix, options, S, &Slen) / BND_DIGIT);
654 case L_BAND_LGAP:
655 return (Int4) ( gband_linear_gap(A + starti - 1, B + startj - 1,
656 endi - starti + 1, endj - startj + 1,
657 matrix, options, S, &Slen) / BND_DIGIT);
658 case L_BAND_QGAP:
659 return (Int4) ( gband_linear_qgap(A + starti - 1, B + startj - 1,
660 endi - starti + 1, endj - startj + 1,
661 matrix, options, S, &Slen) / BND_DIGIT);
662 case L_BAND_L3GAP:
663 return (Int4) ( gband_l3gap(A + starti - 1, B + startj - 1,
664 endi - starti + 1, endj - startj + 1,
665 matrix, options, S, &Slen) / BND_DIGIT);
666 case L_BAND_Q3GAP:
667 return (Int4) ( gband_q3gap(A + starti - 1, B + startj - 1,
668 endi - starti + 1, endj - startj + 1,
669 matrix, options, S, &Slen) / BND_DIGIT);
670 default:
671 ErrPostEx(SEV_ERROR, 0, 0, "Unknown method.");
672 return -1;
673 }
674
675 }
676
677 /**********************************************************************
678 * Gloabal Alignment utility functions
679 **********************************************************************/
680 static Uint1 HSncbi4na_to_blastna[] = {
681 15,/* Gap, 0 */
682 0, /* A, 1 */
683 1, /* C, 2 */
684 6, /* M, 3 */
685 2, /* G, 4 */
686 4, /* R, 5 */
687 9, /* S, 6 */
688 13, /* V, 7 */
689 3, /* T, 8 */
690 8, /* W, 9 */
691 5, /* Y, 10 */
692 12, /* H, 11 */
693 7, /* K, 12 */
694 11, /* D, 13 */
695 10, /* B, 14 */
696 14 /* N, 15 */
697 };
698
699 static Uint1 HSblastna_to_ncbi4na[] = { 1, /* A, 0 */
700 2, /* C, 1 */
701 4, /* G, 2 */
702 8, /* T, 3 */
703 5, /* R, 4 */
704 10, /* Y, 5 */
705 3, /* M, 6 */
706 12, /* K, 7 */
707 9, /* W, 8 */
708 6, /* S, 9 */
709 14, /* B, 10 */
710 13, /* D, 11 */
711 11, /* H, 12 */
712 7, /* V, 13 */
713 15, /* N, 14 */
714 0, /* Gap, 15 */
715 };
716
SetGlobaltOptions(GlobalBandStructPtr gbsp,Int4 lg1_ext,Int4 rg1_ext,Int4 lg2_ext,Int4 rg2_ext,Int4 lg1_open,Int4 lg2_open,Int4 rg1_open,Int4 rg2_open,Int2 gopen,Int2 gext)717 NLM_EXTERN void SetGlobaltOptions(GlobalBandStructPtr gbsp, Int4 lg1_ext, Int4 rg1_ext, Int4 lg2_ext, Int4 rg2_ext, Int4 lg1_open, Int4 lg2_open, Int4 rg1_open, Int4 rg2_open, Int2 gopen, Int2 gext)
718 {
719 PSUGapOptionsPtr opt;
720
721 opt = gbsp->options;
722 opt->gopen = gopen;
723 opt->gext = gext;
724
725 opt->lg1_ext = lg1_ext;
726 opt->rg1_ext = rg1_ext;
727 opt->lg2_ext = lg2_ext;
728 opt->rg2_ext = rg2_ext;
729 opt->lg1_open = lg1_open;
730 opt->lg2_open = lg2_open;
731 opt->rg1_open = rg1_open;
732 opt->rg2_open = rg2_open;
733 }
734
BlastStyleMatCreate(void)735 static Int4Ptr PNTR BlastStyleMatCreate(void){
736 /* Allocate Matrix as one contigous memory block to increase cache-hits&speed
737 */
738 Int4Ptr PNTR matrix;
739 Int4 i;
740 matrix= (Int4Ptr PNTR) Calloc(BLAST_MATRIX_SIZE*BLAST_MATRIX_SIZE,MAX(sizeof(Int4),sizeof(Int4Ptr)));
741 for(i=0;i<BLAST_MATRIX_SIZE-1;i++) {
742 *(matrix+i)=(Int4Ptr) (matrix+BLAST_MATRIX_SIZE+i*BLAST_MATRIX_SIZE);
743 }
744 return matrix;
745 }
746
ChangeGlobalBandMatrix(GlobalBandStructPtr gbsp,Boolean is_prot,CharPtr matrix_name,Int4 penalty,Int4 reward)747 NLM_EXTERN Int2 ChangeGlobalBandMatrix(GlobalBandStructPtr gbsp, Boolean is_prot, CharPtr matrix_name, Int4 penalty, Int4 reward)
748 {
749 BLAST_ScoreBlkPtr sbp;
750 Int4 i, j;
751 Int2 status=0;
752
753 if (is_prot) {
754 sbp = BLAST_ScoreBlkNew((Uint1)Seq_code_ncbistdaa, (Int2)1);
755 sbp->read_in_matrix = TRUE;
756 status = BlastScoreBlkMatFill(sbp, matrix_name);
757 } else {
758 sbp = BLAST_ScoreBlkNew((Uint1)BLASTNA_SEQ_CODE, (Int2)2);
759 sbp->read_in_matrix = FALSE;
760 sbp->penalty = penalty;
761 sbp->reward = reward;
762 status = BlastScoreBlkMatFill(sbp, NULL);
763 }
764 if (status != 0) {
765 ErrPostEx(SEV_WARNING, 0, 0, "BlastScoreBlkMatFill returned non-zero status ");
766 return -1;
767 }
768 if(gbsp->matrix==NULL){
769 gbsp->matrix=(Int4Ptr PNTR) MemFree(gbsp->matrix);
770 }
771 if((gbsp->matrix=BlastStyleMatCreate())==NULL) {
772 ErrPostEx(SEV_WARNING, 0, 0, "BlastStyleMatCreate returned non-zero status ");
773 return -1;
774 }
775
776 if(is_prot) {
777 for (i=0;i<sbp->alphabet_size;i++) {
778 for (j=0;j<sbp->alphabet_size;j++) {
779 gbsp->matrix[i][j] = (Int4) (sbp->matrix[i][j] *BND_DIGIT*(1.000001) );
780 }
781 }
782 } else {
783 for (i=0;i<16;i++) {
784 for (j=0;j<16;j++) {
785 gbsp->matrix[i][j] = (Int4) (sbp->matrix[HSncbi4na_to_blastna[i]][HSncbi4na_to_blastna[j]] *BND_DIGIT*(1.000001) );
786 }
787 }
788 }
789 sbp=BLAST_ScoreBlkDestruct(sbp);
790 return 0;
791 /*
792 BLAST_ScoreBlkPtr sbp;
793 CharPtr name;
794 Int2 i, j, status = 0;
795 Int4Ptr PNTR matrix;
796
797 if (gbsp == NULL) {
798 return 4;
799 }
800 name = (matrix_name) ? matrix_name : "BLOSUM62";
801 if(gbsp->matrix != NULL){
802 gbsp->matrix = (Int4Ptr PNTR ) MemFree(gbsp->matrix);
803 }
804 if (is_prot) {
805 sbp = BLAST_ScoreBlkNew(Seq_code_ncbistdaa, 1);
806 sbp->read_in_matrix = TRUE;
807 BlastScoreSetAmbigRes(sbp, 'X');
808 if((status=BlastScoreBlkMatFill(sbp, name)) != 0) {
809 if((matrix=BlastStyleMatCreate()) == NULL) {
810 sbp=BLAST_ScoreBlkDestruct(sbp);
811 return 0;
812 }
813 for (i=0;i<sbp->alphabet_size;i++) {
814 for (j=0;j<sbp->alphabet_size;j++) {
815 matrix[i][j] = (Int4) (sbp->matrix[i][j] *BND_DIGIT*(1.000001) );
816 }
817 }
818 gbsp->matrix = matrix;
819 }
820 sbp=BLAST_ScoreBlkDestruct(sbp);
821 return status;
822 }
823 sbp = BLAST_ScoreBlkNew(BLASTNA_SEQ_CODE, 2);
824 sbp->read_in_matrix = FALSE;
825 BlastScoreSetAmbigRes(sbp, 'N');
826 sbp->penalty = (Int4) (BND_DIGIT * penalty);
827 sbp->reward = (Int4) (BND_DIGIT * reward);
828 if((status=BlastScoreBlkMatFill(sbp, NULL)) != 0) {
829 if((matrix=BlastStyleMatCreate()) == NULL) {
830 sbp=BLAST_ScoreBlkDestruct(sbp);
831 return 0;
832 }
833 for (i=0;i<16;i++) {
834 for (j=0;j<16;j++) {
835 matrix[i][j] = (Int4) (sbp->matrix[HSncbi4na_to_blastna[i]][HSncbi4na_to_blastna[j]] *BND_DIGIT*(1.000001) );
836 }
837 }
838 gbsp->matrix = matrix;
839 }
840 sbp=BLAST_ScoreBlkDestruct(sbp);
841 return status;
842 */
843 }
844
load_data(SeqLocPtr slp,Boolean is_prot)845 static Uint1Ptr load_data (SeqLocPtr slp,Boolean is_prot)
846 {
847 SeqPortPtr spp;
848 Uint1Ptr seq = NULL;
849 Int4 len;
850 Int4 index;
851 Uint1 code;
852
853 if (is_prot)
854 code = Seq_code_ncbistdaa;
855 else code = Seq_code_ncbi4na;
856
857 len = SeqLocLen(slp);
858 seq = (Uint1Ptr) MemNew(((len)+2)*sizeof(Uint1));
859 spp = SeqPortNewByLoc(slp, code);
860 for (index=0; index < len; index++)
861 {
862 seq[index] = SeqPortGetResidue(spp);
863 }
864 SeqPortFree (spp);
865 return seq;
866 }
867
CreatBandStruct(SeqLocPtr slp1,SeqLocPtr slp2,Int4Ptr PNTR W,Boolean is_prot,Int2 method)868 NLM_EXTERN GlobalBandStructPtr CreatBandStruct(SeqLocPtr slp1, SeqLocPtr slp2, Int4Ptr PNTR W, Boolean is_prot, Int2 method)
869 {
870 GlobalBandStructPtr gbsp;
871 Uint1Ptr seq1,seq2;
872 Int4 len1, len2;
873 Int2 gopen, gext; /* default values */
874 Int4 ma=1, ms = -2; /* default values */
875
876 if (is_prot) {
877 gopen = 5;
878 gext=1;
879 } else {
880 gopen = 10;
881 gext=2;
882 }
883 len1 = SeqLocLen(slp1);
884 len2 = SeqLocLen(slp2);
885
886 seq1 = load_data(slp1, is_prot);
887 seq2 = load_data(slp2, is_prot);
888 gbsp = GlobalBandStructCreate((Uint1)method);
889 gbsp->seq1 = seq1;
890 gbsp->seq2 = seq2;
891 gbsp->seq1_length = len1;
892 gbsp->seq2_length = len2;
893 gbsp->seqloc1 = slp1;
894 gbsp->seqloc2 = slp2;
895
896 if (W != NULL) {
897 gbsp->matrix = W;
898 } else if (is_prot) {
899 ChangeGlobalBandMatrix(gbsp, TRUE, "BLOSUM62", 0, 0);
900 } else {
901 ChangeGlobalBandMatrix(gbsp, FALSE, NULL, ma, ms);
902 }
903
904 /* all end gap penalties are set to zero */
905 SetGlobaltOptions(gbsp, 0, 0, 0, 0, 0, 0, 0, 0, gopen, gext);
906
907 return gbsp;
908 }
909
910 /*
911 Subroutine to define the band to includes segments in the seqalign(s)
912 obtained from Blast. The Band is to be visualized in the space of sequence
913 two (vertical (Y) axis) versus sequence 1 (horizontal (X) axis).
914 By default sequences are filtered for low-complexity sequences.
915 If no alignment is returned, this routine tries again but without filtering.
916
917 the "type" parameter defines how to choose the seqalign(s) to use.
918 type 1 and 2 are meant to give narrow optimal bands.
919 (e.g. when one has to align sequences with an overlap )
920 type = 2 : Only Use the longest Alignment to define the band.
921 type = 1 : Only Use the highest scoring alignment to define the band.
922
923 type = 0: Make the widest possible band that will include all possible
924 blocks found by BLAST. This is meant to globally align
925 two sequences.
926 Additionally, the earliest segment's Y coordinate may increase
927 the upper limit of the band and the last segments may decrease
928 the lower limit of the band. These last two modifications insure
929 that all possible global alignments including those segments are
930 considered, including those with large start or end gaps.
931 */
932
933
SetLowUpFromBlast(PSUGapOptionsPtr opt,Boolean is_prot,Int2 type,Int2 width,SeqLocPtr slp1,SeqLocPtr slp2)934 NLM_EXTERN void SetLowUpFromBlast(PSUGapOptionsPtr opt, Boolean is_prot, Int2 type, Int2 width, SeqLocPtr slp1, SeqLocPtr slp2)
935 {
936 BLAST_OptionsBlkPtr options;
937 SeqAlignPtr seqalign, sap;
938 Int2 index1, i2;
939 Int4 x,y,maxy,miny,mindi,maxdi,mindi_low,maxdi_up;
940 Int4 len,len1,len2,maxlen=0,alnlen,lastx,lowestx,diag;
941 DenseDiagPtr ddp;
942 ScorePtr score;
943 FloatHi val, maxsc = 0;
944 DenseSegPtr dsp;
945 ObjectIdPtr id;
946
947 if (opt == NULL)
948 return;
949 if (width ==0) {
950 width = 30;
951 } else if (width<0) {
952 width=0;
953 }
954 options = BLASTOptionNew((is_prot == FALSE) ? "blastn":"blastp", TRUE);
955 while ((seqalign = BlastTwoSequencesByLoc(slp1, slp2,NULL, options))==NULL && (is_prot==FALSE && options->wordsize>8)) {
956 options->wordsize--;
957 }
958 options=BLASTOptionDelete(options);
959
960 len1 = SeqLocLen(slp1);
961 len2 = SeqLocLen(slp2);
962
963 if (seqalign == NULL || seqalign->segtype > 2) {
964 opt->low = -len2;
965 opt->up = len1;
966 opt->start_diag = opt->low;
967 opt->width = opt->up-opt->low+ 1;
968 return;
969 }
970 miny = len2;
971 maxy = - len1;
972 lastx = 0;
973 lowestx = len1+1;
974 sap = seqalign;
975 for (; sap != NULL; sap = sap->next) {
976 for (score = sap->score; score; score = score->next) {
977 id = score->id;
978 if (StringCmp(id->str, "score") == 0) {
979 if (score->choice == 1) {
980 val = score->value.intvalue;
981 break;
982 }
983 }
984 }
985 if (sap->segtype == 1) {
986 ddp = sap->segs;
987 x = ddp->starts[0];
988 y = ddp->starts[1];
989 alnlen = len = ddp->len;
990 maxdi = mindi = y - x;
991 if(type == 0 ) {
992 if(x+len > lastx) {
993 mindi_low = y-len2+len;
994 lastx = x+len;
995 }
996 if (x<lowestx) {
997 maxdi_up = y;
998 lowestx = x;
999 }
1000 }
1001 } else if (sap->segtype == 2) {
1002 dsp = sap->segs;
1003 x = dsp->starts[0];
1004 y = dsp->starts[1];
1005 alnlen = len = dsp->lens[0];
1006 maxdi = mindi = y - x;
1007 if(type == 0 ) {
1008 if(x+len >= lastx) {
1009 mindi_low = y-len2+len;
1010 lastx = x+len;
1011 }
1012 if (x<lowestx) {
1013 maxdi_up = y;
1014 lowestx = x;
1015 }
1016 }
1017 for (index1=1; index1 < dsp->numseg; index1++) {
1018 i2 = 2 * index1;
1019 x = dsp->starts[i2];
1020 y = dsp->starts[i2+1];
1021 len = dsp->lens[index1];
1022 if (x == -1 || y == -1) {
1023 continue;
1024 }
1025 if (dsp->strands[i2] == Seq_strand_minus ||
1026 dsp->strands[i2+1] == Seq_strand_minus){
1027 continue;
1028 }
1029 diag = y-x;
1030 mindi = MIN(mindi,diag);
1031 maxdi = MAX(maxdi,diag);
1032 alnlen += len;
1033 if(type == 0 ) {
1034 if(x+len > lastx) {
1035 mindi_low = y-len2+len;
1036 lastx = x+len;
1037 }
1038 if (x<lowestx) {
1039 maxdi_up = y;
1040 lowestx = x;
1041 }
1042 }
1043 }
1044 }
1045 switch(type) {
1046 case 0:
1047 mindi-=width;
1048 maxdi+=width;
1049 miny = MIN(miny, mindi);
1050 maxy = MAX(maxy, maxdi);
1051 break;
1052 case 1:
1053 if (maxsc < val) {
1054 maxsc = val;
1055 maxy = maxdi + width;
1056 miny = mindi - width;
1057 }
1058 break;
1059 case 2:
1060 if (maxlen < alnlen) {
1061 maxlen = alnlen;
1062 maxy = maxdi + width;
1063 miny = mindi - width;
1064 }
1065 break;
1066 }
1067 }
1068 if(type==0) {
1069 miny = MIN(miny,mindi_low);
1070 maxy = MAX(maxy,maxdi_up);
1071 }
1072 if (miny<-len1) miny=-len1;
1073 if (maxy>len2) maxy=len2;
1074 opt->low = miny;
1075 opt->up = maxy;
1076 opt->start_diag = miny;
1077 opt->width = maxy - miny + 1;
1078 opt->start_diag=opt->low;
1079 opt->width = opt->up - opt->low + 1;
1080 SeqAlignFree(seqalign);
1081 }
1082
1083
GlobalBandByLoc(GlobalBandStructPtr gbsp,SeqLocPtr slp1,SeqLocPtr slp2,Boolean is_prot,Int2 band_method)1084 NLM_EXTERN SeqAlignPtr GlobalBandByLoc(GlobalBandStructPtr gbsp, SeqLocPtr slp1, SeqLocPtr slp2, Boolean is_prot, Int2 band_method)
1085 {
1086 PSUGapOptionsPtr opt;
1087 SeqAlignPtr seqalign;
1088
1089 opt = gbsp->options;
1090 SetLowUpFromBlast(opt, is_prot, band_method, 0, slp1, slp2);
1091 opt->start_diag = opt->low;
1092 opt->width = opt->up - opt->low + 1;
1093
1094 seqalign = GlobalBandToSeqAlign(gbsp);
1095 AdjustOffSetsInSeqAlign(seqalign, slp1, slp2);
1096 return seqalign;
1097 }
1098
ExtendSeqAlign(SeqAlignPtr sap,Int4 start1,Int4 start2,Int4 stop1,Int4 stop2,Int4 x1,Int4 y1,Int4 x2,Int4 y2)1099 NLM_EXTERN SeqAlignPtr ExtendSeqAlign(SeqAlignPtr sap, Int4 start1, Int4 start2, Int4 stop1, Int4 stop2, Int4 x1, Int4 y1, Int4 x2, Int4 y2)
1100 {
1101 DenseSegPtr dsp;
1102 Int4 index1;
1103 Int4Ptr n_starts, n_lens;
1104 Uint1Ptr n_strands;
1105 Int2 n;
1106 Int2 offset=0;
1107
1108 if (sap==NULL)
1109 return NULL;
1110 if (start1==x1 && start2==y1 && stop1==x2 && stop2==y2)
1111 return sap;
1112 if (sap->segtype == 1) {
1113 DenseDiagPtr ddp;
1114 ddp = sap->segs;
1115 ErrPostEx(SEV_WARNING,0,0,"ExtendSeqAlign doesn't support segtype==1\n");
1116 } else if (sap->segtype == 2) {
1117 dsp = sap->segs;
1118 n = dsp->numseg;
1119 n_starts = (Int4Ptr) MemNew((2*(n+2)+4)*sizeof(Int4));
1120 n_lens = (Int4Ptr) MemNew((n+4)*sizeof(Int4));
1121 n_strands = MemNew(2*(n+2)+4);
1122 if (x1 > start1 || y1 > start2) {
1123 if (x1 > start1) {
1124 n_starts[0] = 0;
1125 n_starts[1] = -1;
1126 n_lens[0] = x1 - start1;
1127 n_strands[0] = Seq_strand_plus;
1128 n_strands[1] = Seq_strand_plus;
1129 } else if (y1 > start2) {
1130 n_starts[0] = -1;
1131 n_starts[1] = 0;
1132 n_lens[0] = y1 - start2;
1133 n_strands[0] = Seq_strand_plus;
1134 n_strands[1] = Seq_strand_plus;
1135 }
1136 offset += 1;
1137 }
1138 for (index1=0; index1<n; index1++) {
1139 n_starts [2*(index1+offset)] = dsp->starts [2*index1];
1140 n_strands[2*(index1+offset)] = dsp->strands[2*index1];
1141 n_starts [2*(index1+offset)+1]=dsp->starts [2*index1+1];
1142 n_strands[2*(index1+offset)+1]=dsp->strands[2*index1+1];
1143 n_lens[index1+offset] = dsp->lens[index1];
1144 }
1145 if (x2 < stop1 || y2 < stop2) {
1146 n += offset;
1147 if (x2 < stop1) {
1148 n_starts[2*n] = x2;
1149 n_starts[2*n+1] = -1;
1150 n_lens[n] = stop1 - x2;
1151 n_strands[2*n] = Seq_strand_unknown;
1152 n_strands[2*n+1] = Seq_strand_plus;
1153 } else if (y2 < stop2) {
1154 n_starts[2*n] = -1;
1155 n_starts[2*n+1] = y2;
1156 n_lens[n] = stop2 - y2;
1157 n_strands[2*n] = Seq_strand_unknown;
1158 n_strands[2*n+1] = Seq_strand_plus;
1159 }
1160 offset += 1;
1161 }
1162 dsp->numseg += offset;
1163 MemFree(dsp->starts);
1164 MemFree(dsp->strands);
1165 MemFree(dsp->lens);
1166 dsp->starts = n_starts;
1167 dsp->lens = n_lens;
1168 dsp->strands = n_strands;
1169 }
1170 return sap;
1171 }
1172
CC_ExtendSeqAlign(SeqAlignPtr sap,Int4 start1,Int4 start2,Int4 stop1,Int4 stop2,Int4 x1,Int4 y1,Int4 x2,Int4 y2,Uint1 strand1,Uint1 strand2)1173 NLM_EXTERN SeqAlignPtr CC_ExtendSeqAlign(SeqAlignPtr sap, Int4 start1, Int4 start2, Int4 stop1, Int4 stop2, Int4 x1, Int4 y1, Int4 x2, Int4 y2, Uint1 strand1, Uint1 strand2)
1174 {
1175 DenseSegPtr dsp;
1176 Int4Ptr n_starts, n_lens;
1177 Uint1Ptr n_strands;
1178 Int4 index1,
1179 minlen = 0;
1180 Int2 n;
1181 Int2 offset=0;
1182
1183 if (sap==NULL)
1184 return NULL;
1185 if (start1==x1 && start2==y1 && stop1==x2 && stop2==y2)
1186 return sap;
1187 if (sap->segtype == 1) {
1188 DenseDiagPtr ddp;
1189 ddp = (DenseDiagPtr) sap->segs;
1190 ErrPostEx(SEV_WARNING,0,0,"ExtendSeqAlign doesn't support segtype==1\n");
1191 } else if (sap->segtype == 2) {
1192 dsp = (DenseSegPtr) sap->segs;
1193 n = dsp->numseg;
1194 n_starts = (Int4Ptr) MemNew((2*(n+2)+4)*sizeof(Int4));
1195 n_lens = (Int4Ptr) MemNew((n+4)*sizeof(Int4));
1196 n_strands = MemNew((2*(n+2)+4)*sizeof(Int1));
1197 if (x1 > start1 && y1 > start2) {
1198 minlen = MIN ((x1-start1), (y1-start2));
1199 n = 0;
1200 }
1201 if ((x1 > start1 || y1 > start2) && ((x1-start1) != (y1-start2))) {
1202 if (x1 > start1) {
1203 n_starts[0] = 0;
1204 n_starts[1] = -1;
1205 n_lens[0] = (x1 - start1) - minlen;
1206 n_strands[0] = strand1;
1207 n_strands[1] = strand2;
1208 }
1209 else if (y1 > start2) {
1210 n_starts[0] = -1;
1211 n_starts[1] = 0;
1212 n_lens[0] = (y1 - start2) - minlen;
1213 n_strands[0] = strand1;
1214 n_strands[1] = strand2;
1215 }
1216 offset += 1;
1217 }
1218 if (minlen > 0) {
1219 n += offset;
1220 n_starts[2*n] = x1 - minlen;
1221 n_starts[2*n+1] = y1 - minlen;
1222 n_lens[n] = minlen;
1223 n_strands[2*n] = strand1;
1224 n_strands[2*n+1] = strand2;
1225 offset += 1;
1226 }
1227 n = dsp->numseg;
1228 for (index1=0; index1<n; index1++) {
1229 n_starts [2*(index1+offset)] = dsp->starts [2*index1];
1230 n_strands[2*(index1+offset)] = dsp->strands[2*index1];
1231 n_starts [2*(index1+offset)+1]=dsp->starts [2*index1+1];
1232 n_strands[2*(index1+offset)+1]=dsp->strands[2*index1+1];
1233 n_lens[index1+offset] = dsp->lens[index1];
1234 }
1235 if (x2 < stop1 && y2 < stop2) {
1236 n += offset;
1237 minlen = MIN ((stop1-x2), (stop2-y2));
1238 n_starts[2*n] = x2 + 1;
1239 n_starts[2*n+1] = y2 +1;
1240 n_lens[n] = minlen;
1241 n_strands[2*n] = strand1;
1242 n_strands[2*n+1] = strand2;
1243 x2 += minlen;
1244 y2 += minlen;
1245 offset += 1;
1246 }
1247 if (x2 < stop1 || y2 < stop2) {
1248 n += offset;
1249 if (x2 < stop1) {
1250 n_starts[2*n] = x2 +1;
1251 n_starts[2*n+1] = -1;
1252 n_lens[n] = stop1 - x2;
1253 n_strands[2*n] = strand1;
1254 n_strands[2*n+1] = strand2;
1255 } else if (y2 < stop2) {
1256 n_starts[2*n] = -1;
1257 n_starts[2*n+1] = y2 +1;
1258 n_lens[n] = stop2 - y2;
1259 n_strands[2*n] = strand1;
1260 n_strands[2*n+1] = strand2;
1261 }
1262 offset += 1;
1263 }
1264 dsp->numseg = n+1;
1265 MemFree(dsp->starts);
1266 MemFree(dsp->strands);
1267 MemFree(dsp->lens);
1268 dsp->starts = n_starts;
1269 dsp->lens = n_lens;
1270 dsp->strands = n_strands;
1271 }
1272 return sap;
1273 }
1274
1275 /* doesn't work for minus strands in SeqAlign */
1276 /* deals with segtypes 1 and 2 only*/
GetAlignExtremes(SeqAlignPtr sap,Int4Ptr xx1,Int4Ptr yy1,Int4Ptr xx2,Int4Ptr yy2)1277 NLM_EXTERN void GetAlignExtremes(SeqAlignPtr sap, Int4Ptr xx1, Int4Ptr yy1, Int4Ptr xx2, Int4Ptr yy2)
1278 {
1279 DenseDiagPtr ddp;
1280 DenseSegPtr dsp;
1281 Int4 index1, x, y, x1, x2, y1, y2;
1282
1283 if (sap == NULL || sap->segtype > 2) {
1284 return;
1285 }
1286 if (sap->segtype == 1) {
1287 ddp = sap->segs;
1288 x1 = ddp->starts[0];
1289 y1 = ddp->starts[1];
1290 x2 = x1 + ddp->len;
1291 y2 = y1 + ddp->len;
1292 } else if (sap->segtype == 2) {
1293 dsp = sap->segs;
1294 x1 = dsp->starts[0];
1295 y1 = dsp->starts[1];
1296 for (index1=0; index1<dsp->numseg; index1++) {
1297 x = dsp->starts[2 * index1];
1298 y = dsp->starts[2 * index1+1];
1299 if (x == -1 || y == -1) {
1300 continue;
1301 }
1302 x2 = x + dsp->lens[index1] - 1;
1303 y2 = y + dsp->lens[index1] - 1;
1304 }
1305 }
1306 *xx1 = x1;
1307 *xx2 = x2;
1308 *yy1 = y1;
1309 *yy2 = y2;
1310 }
1311
1312