← Back to branch summary

~ubuntu-branches/debian/experimental/ncbi-tools6/experimental

« back to all changes in this revision

Viewing changes to algo/blast/api/blast_seqalign.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2005-03-27 12:00:15 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050327120015-embhesp32nj73p9r
Tags: 6.1.20041020-3
http://bugs.debian.org/295110
* Fix FTBFS under GCC 4.0 caused by inconsistent use of "static" on
  functions.  (Closes: #295110.)
* Add a watch file, now that we can.  (Upstream's layout needs version=3.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/api/blast_seqalign.c
 
1
/* $Id: blast_seqalign.c,v 1.45 2004/10/06 15:00:44 dondosha Exp $
 
2
* ===========================================================================
 
3
*
 
4
*                            PUBLIC DOMAIN NOTICE
 
5
*               National Center for Biotechnology Information
 
6
*
 
7
*  This software/database is a "United States Government Work" under the
 
8
*  terms of the United States Copyright Act.  It was written as part of
 
9
*  the author's offical duties as a United States Government employee and
 
10
*  thus cannot be copyrighted.  This software/database is freely available
 
11
*  to the public for use. The National Library of Medicine and the U.S.
 
12
*  Government have not placed any restriction on its use or reproduction.
 
13
*
 
14
*  Although all reasonable efforts have been taken to ensure the accuracy
 
15
*  and reliability of the software and data, the NLM and the U.S.
 
16
*  Government do not and cannot warrant the performance or results that
 
17
*  may be obtained by using this software or data. The NLM and the U.S.
 
18
*  Government disclaim all warranties, express or implied, including
 
19
*  warranties of performance, merchantability or fitness for any particular
 
20
*  purpose.
 
21
*
 
22
*  Please cite the author in any work or product based on this material.
 
23
*
 
24
* ===========================================================================*/
 
25
 
 
26
/*****************************************************************************
 
27
 
 
28
File name: blast_seqalign.c
 
29
 
 
30
Author: Ilya Dondoshansky
 
31
 
 
32
Contents: Conversion of BLAST results to the SeqAlign form
 
33
 
 
34
******************************************************************************
 
35
 * $Revision: 1.45 $
 
36
 * */
 
37
 
 
38
static char const rcsid[] = "$Id: blast_seqalign.c,v 1.45 2004/10/06 15:00:44 dondosha Exp $";
 
39
 
 
40
#include <algo/blast/api/blast_seqalign.h>
 
41
 
 
42
extern Int4 LIBCALL HspArrayPurge PROTO((BlastHSP** hsp_array, 
 
43
                       Int4 hspcnt, Boolean clear_num));
 
44
extern SeqIdPtr GetTheSeqAlignID (SeqIdPtr seq_id);
 
45
extern ScorePtr MakeBlastScore (ScorePtr PNTR old, CharPtr scoretype, 
 
46
                                Nlm_FloatHi prob, Int4 score);
 
47
 
 
48
static ScorePtr 
 
49
GetScoreSetFromBlastHsp(BlastHSP* hsp)
 
50
{
 
51
   ScorePtr     score_set=NULL;
 
52
   double       prob;
 
53
   Int4         score;
 
54
   char*        scoretype;
 
55
 
 
56
   score = hsp->score;
 
57
   if (score > 0)
 
58
      MakeBlastScore(&score_set, "score", 0.0, score);
 
59
 
 
60
   score = hsp->num;
 
61
   scoretype = "sum_n";
 
62
   
 
63
   if (score > 1)
 
64
      MakeBlastScore(&score_set, scoretype, 0.0, score);
 
65
   
 
66
   prob = hsp->evalue;
 
67
   if (hsp->num == 1) {
 
68
      scoretype = "e_value";
 
69
   } else {
 
70
      scoretype = "sum_e";
 
71
   }
 
72
   if (prob >= 0.) {
 
73
      if (prob < 1.0e-180)
 
74
         prob = 0.0;
 
75
      MakeBlastScore(&score_set, scoretype, prob, 0);
 
76
   }
 
77
 
 
78
   /* Calculate bit score from the raw score */
 
79
   if (hsp->bit_score >= 0.)
 
80
      MakeBlastScore(&score_set, "bit_score", hsp->bit_score, 0);
 
81
   
 
82
   if (hsp->num_ident > 0)
 
83
      MakeBlastScore(&score_set, "num_ident", 0.0, hsp->num_ident);
 
84
   
 
85
   return score_set;
 
86
}
 
87
 
 
88
/** Add information on what gis to use to the score set */
 
89
static Int2 AddGiListToScoreSet(ScorePtr score_set, ValNodePtr gi_list)
 
90
 
 
91
{
 
92
   while (gi_list) {
 
93
      MakeBlastScore(&score_set, "use_this_gi", 0.0, gi_list->data.intvalue);
 
94
      gi_list = gi_list->next;
 
95
   }
 
96
   
 
97
   return 0;
 
98
}
 
99
 
 
100
/*************************************************************************
 
101
*
 
102
*       This function fills in the DenseDiag Information from the variable
 
103
*       hsp.  On the first call to this function *old should be
 
104
*       NULL, after that pass in the head of the DenseDiagPtr chain.
 
105
*       The newest DenseDiagPtr is returned.
 
106
*
 
107
************************************************************************/
 
108
 
 
109
static DenseDiagPtr
 
110
BLAST_HSPToDenseDiag(DenseDiagPtr* old, BlastHSP* hsp, Boolean reverse,
 
111
                     Int4 query_length, Int4 subject_length)
 
112
{
 
113
        DenseDiagPtr            ddp, new;
 
114
 
 
115
        new = DenseDiagNew();
 
116
        
 
117
        new->dim = 2;   /* Only 2 is supported in spec. */
 
118
        new->len = hsp->query.length;
 
119
        new->starts = (Int4*) calloc(2, sizeof(Int4));
 
120
        new->strands = (Uint1*) calloc(2, sizeof(Uint1));
 
121
        if (reverse)
 
122
        {
 
123
                if (hsp->subject.frame >= 0)
 
124
                {
 
125
                        new->strands[0] = Seq_strand_plus;
 
126
                        new->starts[0] = hsp->subject.offset;
 
127
                }
 
128
                else
 
129
                {
 
130
                        new->strands[0] = Seq_strand_minus;
 
131
                        new->starts[0] = subject_length - hsp->subject.offset - hsp->subject.length;
 
132
                }
 
133
                if (hsp->query.frame >= 0)
 
134
                {
 
135
                        new->strands[1] = Seq_strand_plus;
 
136
                        new->starts[1] = hsp->query.offset;
 
137
                }
 
138
                else
 
139
                {
 
140
                        new->strands[1] = Seq_strand_minus;
 
141
                        new->starts[1] = query_length - hsp->query.offset - hsp->query.length;
 
142
                }
 
143
        }
 
144
        else
 
145
        {
 
146
                if (hsp->query.frame >= 0)
 
147
                {
 
148
                        new->strands[0] = Seq_strand_plus;
 
149
                        new->starts[0] = hsp->query.offset;
 
150
                }
 
151
                else
 
152
                {
 
153
                        new->strands[0] = Seq_strand_minus;
 
154
                        new->starts[0] = query_length - hsp->query.offset - hsp->query.length;
 
155
                }
 
156
                if (hsp->subject.frame >= 0)
 
157
                {
 
158
                        new->strands[1] = Seq_strand_plus;
 
159
                        new->starts[1] = hsp->subject.offset;
 
160
                }
 
161
                else
 
162
                {
 
163
                        new->strands[1] = Seq_strand_minus;
 
164
                        new->starts[1] = subject_length - hsp->subject.offset - hsp->subject.length;
 
165
                }
 
166
        }
 
167
        new->scores = GetScoreSetFromBlastHsp(hsp);
 
168
 
 
169
/* Go to the end of the chain, and then attach "new" */
 
170
        if (*old)
 
171
        {
 
172
                ddp = *old;
 
173
                while (ddp->next)
 
174
                        ddp = ddp->next;
 
175
                ddp->next = new;
 
176
        }
 
177
        else
 
178
        {
 
179
                *old = new;
 
180
        }
 
181
 
 
182
        new->next = NULL;
 
183
 
 
184
        return new;
 
185
}
 
186
 
 
187
/*************************************************************************
 
188
*
 
189
*       This function fills in the StdSeg Information from the variable
 
190
*       hsp.  On the first call to this function *old should be
 
191
*       NULL, after that pass in the head of the DenseDiagPtr chain.
 
192
*       The newest StdSeg* is returned.
 
193
*
 
194
************************************************************************/
 
195
static StdSeg*
 
196
BLAST_HSPToStdSeg(StdSeg** old, BlastHSP* hsp, Int4 query_length, 
 
197
                  Int4 subject_length, SeqIdPtr sip, Boolean reverse)
 
198
{
 
199
        StdSeg*         ssp,* new;
 
200
        SeqIdPtr                query_sip, subject_sip;
 
201
        SeqIntPtr               seq_int1, seq_int2;
 
202
        SeqLocPtr               slp=NULL;
 
203
 
 
204
        new = StdSegNew();
 
205
/* Duplicate the id and split it up into query and subject parts */
 
206
        query_sip = SeqIdDup(sip);
 
207
        subject_sip = SeqIdDup(sip->next);
 
208
        
 
209
        new->dim = 2;   /* Only 2 is supported in spec. */
 
210
        seq_int1 = SeqIntNew();
 
211
        if (hsp->query.frame == 0)
 
212
        {
 
213
                seq_int1->from = hsp->query.offset;
 
214
                seq_int1->to = hsp->query.offset + hsp->query.length - 1;
 
215
                seq_int1->strand = Seq_strand_unknown;
 
216
        }
 
217
        else if (hsp->query.frame < 0)
 
218
        {
 
219
                seq_int1->to = query_length - CODON_LENGTH*hsp->query.offset
 
220
         + hsp->query.frame;
 
221
                seq_int1->from = query_length
 
222
         - CODON_LENGTH*(hsp->query.offset+hsp->query.length)
 
223
         + hsp->query.frame + 1;
 
224
                seq_int1->strand = Seq_strand_minus;
 
225
        }
 
226
        else if (hsp->query.frame > 0)
 
227
        {
 
228
                seq_int1->from = CODON_LENGTH*(hsp->query.offset) + hsp->query.frame - 1;
 
229
                seq_int1->to = CODON_LENGTH*(hsp->query.offset+hsp->query.length)
 
230
         + hsp->query.frame - 2;
 
231
                seq_int1->strand = Seq_strand_plus;
 
232
        }
 
233
        seq_int1->id = query_sip;
 
234
        seq_int2 = SeqIntNew();
 
235
        if (hsp->subject.frame == 0)
 
236
        {
 
237
                seq_int2->from = hsp->subject.offset;
 
238
                seq_int2->to = hsp->subject.offset + hsp->subject.length - 1;
 
239
                seq_int2->strand = Seq_strand_unknown;
 
240
        } 
 
241
        else if (hsp->subject.frame < 0)
 
242
        {
 
243
                seq_int2->from = subject_length - 
 
244
         CODON_LENGTH*(hsp->subject.offset + hsp->subject.length) 
 
245
         + hsp->subject.frame + 1;
 
246
                seq_int2->to = subject_length - CODON_LENGTH*(hsp->subject.offset) 
 
247
         + hsp->subject.frame;
 
248
                seq_int2->strand = Seq_strand_minus;
 
249
        }
 
250
        else if (hsp->subject.frame > 0)
 
251
        {
 
252
                seq_int2->from = CODON_LENGTH*(hsp->subject.offset) + 
 
253
         hsp->subject.frame - 1;
 
254
                seq_int2->to = CODON_LENGTH*(hsp->subject.offset + hsp->subject.length)
 
255
         + hsp->subject.frame - 2;
 
256
                seq_int2->strand = Seq_strand_plus;
 
257
        }
 
258
        seq_int2->id = subject_sip;
 
259
 
 
260
        if (reverse)
 
261
        {
 
262
                ValNodeAddPointer(&slp, SEQLOC_INT, seq_int2); 
 
263
                ValNodeAddPointer(&slp, SEQLOC_INT, seq_int1); 
 
264
        }
 
265
        else
 
266
        {
 
267
                ValNodeAddPointer(&slp, SEQLOC_INT, seq_int1); 
 
268
                ValNodeAddPointer(&slp, SEQLOC_INT, seq_int2); 
 
269
        }
 
270
        new->loc = slp;
 
271
 
 
272
        new->scores = GetScoreSetFromBlastHsp(hsp);
 
273
 
 
274
/* Go to the end of the chain, and then attach "new" */
 
275
        if (*old)
 
276
        {
 
277
                ssp = *old;
 
278
                while (ssp->next)
 
279
                        ssp = ssp->next;
 
280
                ssp->next = new;
 
281
        }
 
282
        else
 
283
        {
 
284
                *old = new;
 
285
        }
 
286
 
 
287
        new->next = NULL;
 
288
 
 
289
        return new;
 
290
}
 
291
 
 
292
/************************************************************************
 
293
*
 
294
*       This function assembles all the components of the Seq-align from
 
295
*       a "sparse" BLAST HitList.  "sparse" means that the hitlist 
 
296
*       may contain no sequence and not even a descriptor.  It is only 
 
297
*       required to contain the sequence_number that readdb refers to
 
298
*       and scoring/alignment information.
 
299
*
 
300
*       If dbname is non-NULL, then only a general ("gnl") ID is 
 
301
*       issued, with the ordinal number of the subject sequence in
 
302
*       the ObjectIdPtr.
 
303
*
 
304
*       Boolean reverse: reverse the query and db order in SeqAlign.
 
305
*
 
306
************************************************************************/
 
307
static Int2 
 
308
BLAST_UngappedHSPToSeqAlign(EBlastProgramType program_number, 
 
309
   BlastHSPList* hsp_list, SeqIdPtr query_id, 
 
310
   SeqIdPtr subject_id, Int4 query_length,
 
311
   Int4 subject_length, SeqAlignPtr* seqalign_ptr)
 
312
{
 
313
   BlastHSP* hsp;
 
314
   DenseDiagPtr ddp_head=NULL, ddp;
 
315
   SeqIdPtr sip;
 
316
   SeqIdPtr new_sip=NULL;
 
317
   StdSeg* ssp_head=NULL,* ssp;
 
318
   SeqAlignPtr seqalign;
 
319
   Int4 hsp_cnt, index2, hspset_cnt_old;
 
320
   Boolean getdensediag = 
 
321
      (program_number == eBlastTypeBlastn ||
 
322
       program_number == eBlastTypeRpsBlast ||
 
323
       program_number == eBlastTypeBlastp);
 
324
 
 
325
        ddp_head = NULL;
 
326
        ssp_head = NULL;
 
327
        sip = NULL;
 
328
 
 
329
 
 
330
   seqalign = SeqAlignNew();
 
331
   seqalign->type = 2;          /* alignment is diags */
 
332
 
 
333
   hspset_cnt_old = -1;
 
334
   hsp_cnt = hsp_list->hspcnt;
 
335
 
 
336
   for (index2=0; index2<hsp_cnt; index2++) {
 
337
      hsp = hsp_list->hsp_array[index2];
 
338
 
 
339
      sip = GetTheSeqAlignID(query_id);
 
340
      sip->next = SeqIdDup(subject_id);
 
341
 
 
342
      if (getdensediag) {
 
343
         ddp = BLAST_HSPToDenseDiag(&ddp_head, hsp, FALSE, query_length, 
 
344
                                    subject_length);
 
345
         ddp->id = sip;
 
346
      } else {
 
347
         ssp = BLAST_HSPToStdSeg(&ssp_head, hsp, query_length, 
 
348
                                 subject_length, sip, FALSE);
 
349
         ssp->ids = sip;
 
350
      }
 
351
      sip = NULL; /* This SeqIdPtr is now on the SeqAlign. */
 
352
   }
 
353
 
 
354
   if (getdensediag) {
 
355
      seqalign->segs = ddp_head;
 
356
      seqalign->segtype = 1;  /* DenseDiag */
 
357
   } else {
 
358
      seqalign->segs = ssp_head;
 
359
      seqalign->segtype = 3;  /* StdSeg */
 
360
   }
 
361
 
 
362
   if (new_sip)
 
363
      new_sip = SeqIdFree(new_sip);
 
364
 
 
365
   *seqalign_ptr = seqalign;
 
366
 
 
367
   return 0;
 
368
}
 
369
 
 
370
/** Get the current position. */
 
371
static Int4 get_current_pos(Int4* pos, Int4 length)
 
372
{
 
373
    Int4 val;
 
374
    if(*pos < 0)
 
375
        val = -(*pos + length -1);
 
376
    else
 
377
        val = *pos;
 
378
    *pos += length;
 
379
    return val;
 
380
}
 
381
 
 
382
/** Given gap information in an edit block form, fill the starts, lengths and
 
383
 * strands arrays
 
384
 */
 
385
Boolean GapCollectDataForSeqalign(GapEditBlock* edit_block, 
 
386
                                  GapEditScript* curr_in, Int4 numseg,
 
387
                                  Int4** start_out, 
 
388
                                  Int4** length_out,
 
389
                                  Uint1** strands_out,
 
390
                                  Int4* start1, Int4* start2)
 
391
{
 
392
    GapEditScript* curr;
 
393
    Boolean reverse, translate1, translate2;
 
394
    Int2 frame1, frame2;
 
395
    Int4 begin1, begin2, index, length1, length2;
 
396
    Int4 original_length1, original_length2, i;
 
397
    Int4* length,* start;
 
398
    Uint1 strand1, strand2;
 
399
    Uint1* strands;
 
400
    
 
401
    reverse = edit_block->reverse;      
 
402
    length1 = edit_block->length1;
 
403
    length2 = edit_block->length2;
 
404
    original_length1 = edit_block->original_length1;
 
405
    original_length2 = edit_block->original_length2;
 
406
    translate1 = edit_block->translate1;
 
407
    translate2 = edit_block->translate2;
 
408
    frame1 = edit_block->frame1;
 
409
    frame2 = edit_block->frame2;
 
410
    
 
411
    if (frame1 > 0)
 
412
        strand1 = Seq_strand_plus; 
 
413
    else if (frame1 < 0)
 
414
        strand1 = Seq_strand_minus; 
 
415
    else
 
416
        strand1 = Seq_strand_unknown; 
 
417
    
 
418
    if (frame2 > 0)
 
419
        strand2 = Seq_strand_plus; 
 
420
    else if (frame2 < 0)
 
421
        strand2 = Seq_strand_minus; 
 
422
    else
 
423
        strand2 = Seq_strand_unknown; 
 
424
 
 
425
    start = (Int4 *) calloc((2*numseg+1), sizeof(Int4));
 
426
    length = (Int4 *) calloc((numseg+1), sizeof(Int4));
 
427
    strands = (Uint1 *) calloc((2*numseg+1), sizeof(Uint1));
 
428
 
 
429
    index=0;
 
430
    for (i = 0, curr=curr_in; curr && i < numseg; curr=curr->next, i++) {
 
431
        switch(curr->op_type) {
 
432
        case eGapAlignDecline:
 
433
        case eGapAlignSub:
 
434
            if (strand1 != Seq_strand_minus) {
 
435
                if(translate1 == FALSE)
 
436
                    begin1 = get_current_pos(start1, curr->num);
 
437
                else
 
438
                    begin1 = frame1 - 1 + CODON_LENGTH*get_current_pos(start1, curr->num);
 
439
            } else {
 
440
                if(translate1 == FALSE)
 
441
                    begin1 = length1 - get_current_pos(start1, curr->num) - curr->num;
 
442
                else
 
443
                    begin1 = original_length1 - CODON_LENGTH*(get_current_pos(start1, curr->num)+curr->num) + frame1 + 1;
 
444
            }
 
445
                                        
 
446
            if (strand2 != Seq_strand_minus) {
 
447
                if(translate2 == FALSE)
 
448
                    begin2 = get_current_pos(start2, curr->num);
 
449
                else
 
450
                    begin2 = frame2 - 1 + CODON_LENGTH*get_current_pos(start2, curr->num);
 
451
            } else {
 
452
                if(translate2 == FALSE)
 
453
                    begin2 = length2 - get_current_pos(start2, curr->num) - curr->num;
 
454
                else
 
455
                    begin2 = original_length2 - CODON_LENGTH*(get_current_pos(start2, curr->num)+curr->num) + frame2 + 1;
 
456
            }
 
457
            
 
458
            if (reverse) {
 
459
                strands[2*index] = strand2;
 
460
                strands[2*index+1] = strand1;
 
461
                start[2*index] = begin2;
 
462
                start[2*index+1] = begin1;
 
463
            } else {
 
464
                strands[2*index] = strand1;
 
465
                strands[2*index+1] = strand2;
 
466
                start[2*index] = begin1;
 
467
                start[2*index+1] = begin2;
 
468
            }
 
469
            
 
470
            break;
 
471
            
 
472
        case eGapAlignDel:
 
473
            begin1 = -1;
 
474
            if (strand2 != Seq_strand_minus) {
 
475
                if(translate2 == FALSE)
 
476
                    begin2 = get_current_pos(start2, curr->num);
 
477
                else
 
478
                    begin2 = frame2 - 1 + CODON_LENGTH*get_current_pos(start2, curr->num);
 
479
            } else {
 
480
                if(translate2 == FALSE)
 
481
                    begin2 = length2 - get_current_pos(start2, curr->num) - curr->num;
 
482
                else
 
483
                    begin2 = original_length2 - CODON_LENGTH*(get_current_pos(start2, curr->num)+curr->num) + frame2 + 1;
 
484
            }
 
485
            
 
486
            if (reverse) {
 
487
                strands[2*index] = strand2;
 
488
                if (index > 0)
 
489
                    strands[2*index+1] = strands[2*(index-1)+1];
 
490
                else
 
491
                    strands[2*index+1] = Seq_strand_unknown;
 
492
                start[2*index] = begin2;
 
493
                start[2*index+1] = begin1;
 
494
            } else {
 
495
                if (index > 0)
 
496
                    strands[2*index] = strands[2*(index-1)];
 
497
                else
 
498
                    strands[2*index] = Seq_strand_unknown;
 
499
                strands[2*index+1] = strand2;
 
500
                start[2*index] = begin1;
 
501
                start[2*index+1] = begin2;
 
502
            }
 
503
            
 
504
            break;
 
505
            
 
506
        case eGapAlignIns:
 
507
            if (strand1 != Seq_strand_minus) {
 
508
                if(translate1 == FALSE)
 
509
                    begin1 = get_current_pos(start1, curr->num);
 
510
                else
 
511
                    begin1 = frame1 - 1 + CODON_LENGTH*get_current_pos(start1, curr->num);
 
512
            } else {
 
513
                if(translate1 == FALSE)
 
514
                    begin1 = length1 - get_current_pos(start1, curr->num) - curr->num;
 
515
                else
 
516
                    begin1 = original_length1 - CODON_LENGTH*(get_current_pos(start1, curr->num)+curr->num) + frame1 + 1;
 
517
            }
 
518
            begin2 = -1;
 
519
            if (reverse) {
 
520
                if (index > 0)
 
521
                    strands[2*index] = strands[2*(index-1)];
 
522
                else
 
523
                    strands[2*index] = Seq_strand_unknown;
 
524
                strands[2*index+1] = strand1;
 
525
                start[2*index] = begin2;
 
526
                start[2*index+1] = begin1;
 
527
            } else {
 
528
                strands[2*index] = strand1;
 
529
                if (index > 0)
 
530
                    strands[2*index+1] = strands[2*(index-1)+1];
 
531
                else
 
532
                    strands[2*index+1] = Seq_strand_unknown;
 
533
                start[2*index] = begin1;
 
534
                start[2*index+1] = begin2;
 
535
            }
 
536
            
 
537
            break;
 
538
        default:
 
539
            break;
 
540
        }
 
541
        length[index] = curr->num;
 
542
        index++;
 
543
    }    
 
544
 
 
545
    if (start_out)
 
546
       *start_out = start;
 
547
    else
 
548
       sfree(start);
 
549
    if (length_out)
 
550
       *length_out = length;
 
551
    else
 
552
       sfree(length);
 
553
    if (strands_out)
 
554
       *strands_out = strands;
 
555
    else
 
556
       sfree(strands);
 
557
 
 
558
    return TRUE;
 
559
}
 
560
 
 
561
static void GapCorrectUASequence(GapEditBlock* edit_block)
 
562
{
 
563
    GapEditScript* curr,* curr_last,* curr_last2;
 
564
    Boolean last_indel;
 
565
 
 
566
    last_indel = FALSE;
 
567
    curr_last = NULL;
 
568
    curr_last2 = NULL;
 
569
 
 
570
    for (curr=edit_block->esp; curr; curr = curr->next) {
 
571
 
 
572
        if(curr->op_type == eGapAlignDecline && last_indel == TRUE) {
 
573
            /* This is invalid condition and regions should be
 
574
               exchanged */
 
575
 
 
576
            if(curr_last2 != NULL)
 
577
                curr_last2->next = curr;
 
578
            else
 
579
                edit_block->esp = curr; /* First element in a list */
 
580
            
 
581
            curr_last->next = curr->next;
 
582
            curr->next = curr_last;
 
583
        }
 
584
        
 
585
        last_indel = FALSE;
 
586
        
 
587
        if(curr->op_type == eGapAlignIns || curr->op_type == eGapAlignDel) {
 
588
            last_indel = TRUE;
 
589
            curr_last2 = curr_last;
 
590
        }
 
591
 
 
592
        curr_last = curr;
 
593
    }
 
594
    return;
 
595
}
 
596
 
 
597
static SeqAlignPtr GapMakeSeqAlign(SeqIdPtr query_id, SeqIdPtr subject_id, 
 
598
                                   Boolean reverse, Boolean translate1, 
 
599
                                   Boolean translate2, Int4 numseg,
 
600
                                   Int4* length, Int4* start, 
 
601
                                   Uint1* strands)
 
602
{
 
603
    SeqAlignPtr sap;
 
604
    DenseSeg* dsp;
 
605
    StdSeg* sseg,* sseg_head,* sseg_old;
 
606
    SeqLocPtr slp, slp1, slp2;
 
607
    SeqIntPtr seq_int1;
 
608
    Int4 index;
 
609
    Boolean tmp_value;
 
610
 
 
611
    sap = SeqAlignNew();
 
612
    
 
613
    sap->dim =2; /**only two dimention alignment**/
 
614
    
 
615
    /**make the Denseg Object for SeqAlign**/
 
616
    if (translate1 == FALSE && translate2 == FALSE) {
 
617
        sap->segtype = SAS_DENSEG; /** use denseg to store the alignment **/
 
618
        sap->type = SAT_PARTIAL;   /**partial for gapped translating search.*/
 
619
        dsp = DenseSegNew();
 
620
        dsp->dim = 2;
 
621
        dsp->numseg = numseg;
 
622
        if (reverse) {
 
623
            dsp->ids = SeqIdDup(subject_id);
 
624
            dsp->ids->next = SeqIdDup(query_id);
 
625
        } else {
 
626
            dsp->ids = SeqIdDup(query_id);
 
627
            dsp->ids->next = SeqIdDup(subject_id);
 
628
        }
 
629
        dsp->starts = start;
 
630
        dsp->strands = strands;
 
631
        dsp->lens = length;
 
632
        sap->segs = dsp;
 
633
        sap->next = NULL;
 
634
    } else { /****/
 
635
        sap->type = SAT_PARTIAL; /**partial for gapped translating search. */
 
636
        sap->segtype = SAS_STD;  /**use stdseg to store the alignment**/
 
637
        sseg_head = NULL;
 
638
        sseg_old = NULL;
 
639
 
 
640
        if(reverse) {  /* Reverting translation flags */
 
641
            tmp_value = translate1;
 
642
            translate1 = translate2;
 
643
            translate2 = tmp_value;
 
644
        }
 
645
        
 
646
        for (index=0; index<numseg; index++) {
 
647
            sseg = StdSegNew();
 
648
            sseg->dim = 2;
 
649
            if (sseg_head == NULL) {
 
650
                sseg_head = sseg;
 
651
            }
 
652
            if (reverse) {
 
653
                sseg->ids = SeqIdDup(subject_id);
 
654
                sseg->ids->next = SeqIdDup(query_id);
 
655
            } else {
 
656
                sseg->ids = SeqIdDup(query_id);
 
657
                sseg->ids->next = SeqIdDup(subject_id);
 
658
            }
 
659
 
 
660
            slp1 = NULL;
 
661
            if (start[2*index] != -1) {
 
662
                seq_int1 = SeqIntNew();
 
663
                seq_int1->from = start[2*index];
 
664
                if (translate1)
 
665
                    seq_int1->to = start[2*index] + CODON_LENGTH*length[index] - 1;
 
666
                else
 
667
                    seq_int1->to = start[2*index] + length[index] - 1;
 
668
                seq_int1->strand = strands[2*index];
 
669
 
 
670
                if(reverse) { 
 
671
                    seq_int1->id = SeqIdDup(subject_id);
 
672
                } else {
 
673
                    seq_int1->id = SeqIdDup(query_id);
 
674
                }
 
675
 
 
676
                ValNodeAddPointer(&slp1, SEQLOC_INT, seq_int1);
 
677
            } else {
 
678
                if(reverse) { 
 
679
                    ValNodeAddPointer(&slp1, SEQLOC_EMPTY, SeqIdDup(subject_id));
 
680
                } else {
 
681
                    ValNodeAddPointer(&slp1, SEQLOC_EMPTY, SeqIdDup(query_id));
 
682
                }
 
683
            }
 
684
            slp2 = NULL;
 
685
            if (start[2*index+1] != -1) {
 
686
                seq_int1 = SeqIntNew();
 
687
                seq_int1->from = start[2*index+1];
 
688
                if (translate2)
 
689
                    seq_int1->to = start[2*index+1] + CODON_LENGTH*length[index] - 1;
 
690
                else
 
691
                    seq_int1->to = start[2*index+1] + length[index] - 1;
 
692
                seq_int1->strand = strands[2*index+1];
 
693
 
 
694
                if(reverse) { 
 
695
                    seq_int1->id = SeqIdDup(query_id);
 
696
                } else {
 
697
                    seq_int1->id = SeqIdDup(subject_id);
 
698
                }
 
699
 
 
700
                ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int1);
 
701
            } else {
 
702
                if(reverse) { 
 
703
                    ValNodeAddPointer(&slp2, SEQLOC_EMPTY, SeqIdDup(query_id));
 
704
                } else {
 
705
                    ValNodeAddPointer(&slp2, SEQLOC_EMPTY, SeqIdDup(subject_id));
 
706
                }
 
707
            }
 
708
            /*
 
709
              if (reverse) {
 
710
              slp = slp2;
 
711
              slp2->next = slp1;
 
712
              } else {
 
713
              slp = slp1;
 
714
              slp1->next = slp2;
 
715
              } */
 
716
 
 
717
            slp = slp1;
 
718
            slp1->next = slp2;
 
719
 
 
720
            sseg->loc = slp;
 
721
            
 
722
            if (sseg_old)
 
723
                sseg_old->next = sseg;
 
724
            sseg_old = sseg;
 
725
        }
 
726
        sap->segs = sseg_head;
 
727
        sap->next = NULL;
 
728
        
 
729
        sfree(start);
 
730
        sfree(length);
 
731
        sfree(strands);
 
732
    }
 
733
 
 
734
    return sap;
 
735
}
 
736
 
 
737
/** Convert an EditScript chain to a SeqAlign of type DenseSeg.
 
738
 * Used for a non-simple interval (i.e., one without subs. or 
 
739
 * deletions).  
 
740
 * 
 
741
 * The first SeqIdPtr in the chain of subject_id and query_id is duplicated for
 
742
 * the SeqAlign.
 
743
 */
 
744
SeqAlignPtr LIBCALL
 
745
GapEditBlockToSeqAlign(GapEditBlock* edit_block, SeqIdPtr subject_id, SeqIdPtr query_id)
 
746
 
 
747
{
 
748
    GapEditScript* curr,* curr_last;
 
749
    Int4 numseg, start1, start2;
 
750
    Int4* length,* start;
 
751
    Uint1* strands;
 
752
    Boolean is_disc_align = FALSE;
 
753
    SeqAlignPtr sap, sap_disc, sap_head, sap_tail;
 
754
    Boolean skip_region = FALSE;
 
755
 
 
756
    is_disc_align = FALSE; numseg = 0;
 
757
 
 
758
    for (curr=edit_block->esp; curr; curr = curr->next) {
 
759
        numseg++;
 
760
        if(/*edit_block->discontinuous && */curr->op_type == eGapAlignDecline)
 
761
           is_disc_align = TRUE;
 
762
    }
 
763
    
 
764
    start1 = edit_block->start1;
 
765
    start2 = edit_block->start2;
 
766
    
 
767
    /* If no eGapAlignDecline regions exists output seqalign will be
 
768
       regular Den-Seg or Std-seg */
 
769
    if(is_disc_align == FALSE) {
 
770
        /* Please note, that edit_block passed only for data like
 
771
           strand, translate, reverse etc. Real data is taken starting
 
772
           from "curr" and taken only "numseg" segments */
 
773
        
 
774
        GapCollectDataForSeqalign(edit_block, edit_block->esp, numseg,
 
775
                                  &start, &length, &strands, &start1, &start2);
 
776
        
 
777
        /* Result of this function will be either den-seg or Std-seg
 
778
           depending on translation options */
 
779
        sap = GapMakeSeqAlign(query_id, subject_id, edit_block->reverse, 
 
780
                              edit_block->translate1, edit_block->translate2, 
 
781
                              numseg, length, start, strands);
 
782
    } else {
 
783
 
 
784
        /* By request of Steven Altschul - we need to have 
 
785
           the unaligned part being to the left if it is adjacent to the
 
786
           gap (insertion or deletion) - so this function will do
 
787
           shaffeling */
 
788
 
 
789
        GapCorrectUASequence(edit_block); 
 
790
 
 
791
        sap_disc = SeqAlignNew();
 
792
        sap_disc->dim = 2;
 
793
        sap_disc->type = SAT_PARTIAL; /* ordered segments, over part of seq */
 
794
        sap_disc->segtype = SAS_DISC; /* discontinuous alignment */
 
795
        
 
796
        curr_last = edit_block->esp; 
 
797
        sap_head = NULL; sap_tail = NULL;
 
798
        while(curr_last != NULL) {
 
799
            skip_region = FALSE;                        
 
800
            for (numseg = 0, curr = curr_last; 
 
801
                 curr; curr = curr->next, numseg++) {
 
802
 
 
803
                if(curr->op_type == eGapAlignDecline) {
 
804
                    if(numseg != 0) { /* End of aligned area */
 
805
                        break;
 
806
                    } else {
 
807
                        while(curr && curr->op_type == eGapAlignDecline) {
 
808
                            numseg++;
 
809
                            curr = curr->next; 
 
810
                        }
 
811
                        skip_region = TRUE;                        
 
812
                        break;
 
813
                    }
 
814
                }
 
815
            }
 
816
 
 
817
            GapCollectDataForSeqalign(edit_block, curr_last, numseg,
 
818
                                      &start, &length, &strands, 
 
819
                                      &start1, &start2);
 
820
            
 
821
            if(!skip_region) {            
 
822
                sap = GapMakeSeqAlign(query_id, subject_id, 
 
823
                                      edit_block->reverse, 
 
824
                                      edit_block->translate1, 
 
825
                                      edit_block->translate2,
 
826
                                      numseg, length, start, strands);
 
827
                
 
828
                /* Collecting all seqaligns into single linked list */
 
829
                if(sap_tail == NULL) {
 
830
                    sap_head = sap_tail = sap;
 
831
                } else {
 
832
                    sap_tail->next = sap;
 
833
                    sap_tail = sap;
 
834
                }
 
835
            }
 
836
            
 
837
            curr_last = curr;
 
838
        }
 
839
        sap_disc->segs = sap_head;
 
840
        sap = sap_disc;
 
841
    }
 
842
 
 
843
    return sap;
 
844
}
 
845
 
 
846
/** This function is used for Out-Of-Frame traceback conversion
 
847
 * Convert an OOF EditScript chain to a SeqAlign of type StdSeg.
 
848
 * Used for a non-simple interval (i.e., one without subs. or 
 
849
 * deletions).  
 
850
 * The first SeqIdPtr in the chain of subject_id and query_id is 
 
851
 * duplicated for the SeqAlign.
 
852
*/
 
853
static SeqAlignPtr LIBCALL
 
854
BLAST_OOFEditBlockToSeqAlign(EBlastProgramType program, 
 
855
   GapEditBlock* edit_block, SeqIdPtr query_id, SeqIdPtr subject_id)
 
856
{
 
857
    Boolean reverse = FALSE;
 
858
    GapEditScript* curr,* esp;
 
859
    Int2 frame1, frame2;
 
860
    Int4 start1, start2;
 
861
    Int4 original_length1, original_length2;
 
862
    SeqAlignPtr sap;
 
863
    SeqIntPtr seq_int1, seq_int2;
 
864
    SeqIntPtr seq_int1_last = NULL, seq_int2_last = NULL;
 
865
    SeqIdPtr sip, id1, id2;
 
866
    SeqLocPtr slp, slp1, slp2;
 
867
    StdSeg* sseg,* sseg_head,* sseg_old;
 
868
    Uint1 strand1, strand2;
 
869
    Boolean first_shift;
 
870
 
 
871
    if (program == eBlastTypeBlastx) {
 
872
       reverse = TRUE;
 
873
       start1 = edit_block->start2;
 
874
       start2 = edit_block->start1;
 
875
       frame1 = edit_block->frame2;
 
876
       frame2 = edit_block->frame1;
 
877
       original_length1 = edit_block->original_length2;
 
878
       original_length2 = edit_block->original_length1;
 
879
       id1 = subject_id;
 
880
       id2 = query_id;
 
881
    } else { 
 
882
       start1 = edit_block->start1;
 
883
       start2 = edit_block->start2;
 
884
       frame1 = edit_block->frame1;
 
885
       frame2 = edit_block->frame2;
 
886
       original_length1 = edit_block->original_length1;
 
887
       original_length2 = edit_block->original_length2;
 
888
       id1 = query_id;
 
889
       id2 = subject_id;
 
890
    }
 
891
 
 
892
    if(frame1 > 0) 
 
893
        strand1 = Seq_strand_plus;
 
894
    else if (frame1 < 0)
 
895
        strand1 = Seq_strand_minus;
 
896
    else
 
897
        strand1 = Seq_strand_unknown;
 
898
    
 
899
    if(frame2 > 0) 
 
900
        strand2 = Seq_strand_plus;
 
901
    else if (frame2 < 0)
 
902
        strand2 = Seq_strand_minus;
 
903
    else
 
904
        strand2 = Seq_strand_unknown;
 
905
    
 
906
    esp = edit_block->esp;
 
907
    
 
908
    sap = SeqAlignNew();
 
909
    
 
910
    sap->dim =2; /**only two dimention alignment**/
 
911
    
 
912
    sap->type =3; /**partial for gapped translating search. */
 
913
    sap->segtype =3; /**use denseg to store the alignment**/
 
914
    sseg_head = NULL;
 
915
    sseg_old = NULL;
 
916
    esp = edit_block->esp;
 
917
 
 
918
    first_shift = FALSE;
 
919
 
 
920
    for (curr=esp; curr; curr=curr->next) {
 
921
 
 
922
        slp1 = NULL;
 
923
        slp2 = NULL;
 
924
        
 
925
        switch (curr->op_type) {
 
926
        case eGapAlignDel: /* deletion of three nucleotides. */
 
927
            
 
928
            first_shift = FALSE;
 
929
 
 
930
            seq_int1 = SeqIntNew();
 
931
            seq_int1->from = get_current_pos(&start1, curr->num);
 
932
            seq_int1->to = start1 - 1;            
 
933
 
 
934
            if(seq_int1->to >= original_length1)
 
935
                seq_int1->to = original_length1-1;
 
936
            
 
937
            seq_int1->id = SeqIdDup(id1);
 
938
            seq_int1->strand = strand1;
 
939
 
 
940
            ValNodeAddPointer(&slp1, SEQLOC_INT, seq_int1);
 
941
 
 
942
            /* Empty nucleotide piece */
 
943
            ValNodeAddPointer(&slp2, SEQLOC_EMPTY, SeqIdDup(id2));
 
944
            
 
945
            seq_int1_last = seq_int1;
 
946
            /* Keep previous seq_int2_last, in case there is a frame shift
 
947
               immediately after this gap */
 
948
            
 
949
            break;
 
950
 
 
951
        case eGapAlignIns: /* insertion of three nucleotides. */
 
952
 
 
953
            /* If gap is followed after frameshift - we have to
 
954
               add this element for the alignment to be correct */
 
955
            
 
956
            if(first_shift == TRUE) { /* Second frameshift in a row */
 
957
                /* Protein coordinates */
 
958
                seq_int1 = SeqIntNew();
 
959
                seq_int1->from =  get_current_pos(&start1, 1);
 
960
                seq_int1->to = start1 - 1;
 
961
 
 
962
                if(seq_int1->to >= original_length1)
 
963
                    seq_int1->to = original_length1-1;
 
964
                
 
965
                seq_int1->id = SeqIdDup(id1);
 
966
                seq_int1->strand = strand1;
 
967
                
 
968
                ValNodeAddPointer(&slp1, SEQLOC_INT, seq_int1);
 
969
                
 
970
                /* Nucleotide scale shifted by op_type */
 
971
                seq_int2 = SeqIntNew();
 
972
 
 
973
                seq_int2->from = get_current_pos(&start2, 3);
 
974
                seq_int2->to = start2 - 1;
 
975
 
 
976
                if(seq_int2->to >= original_length2) {
 
977
                    seq_int2->to = original_length2 - 1;
 
978
                    seq_int1->to--;
 
979
                }
 
980
 
 
981
                /* Transfer to DNA minus strand coordinates */
 
982
                if(strand2 == Seq_strand_minus) {
 
983
                    int tmp_int;
 
984
                    tmp_int = seq_int2->to;
 
985
                    seq_int2->to = original_length2 - seq_int2->from - 1;
 
986
                    seq_int2->from = original_length2 - tmp_int - 1;
 
987
                }
 
988
            
 
989
                seq_int2->id = SeqIdDup(id2);
 
990
                seq_int2->strand = strand2;
 
991
                
 
992
                ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int2);
 
993
 
 
994
                /* seq_int1_last = seq_int1; 
 
995
                   seq_int2_last = seq_int2; */
 
996
 
 
997
                /* first_shift = FALSE; */
 
998
 
 
999
                if (reverse) {
 
1000
                    slp = slp2;
 
1001
                    slp2->next = slp1;
 
1002
                    sip = SeqIdDup(id2);
 
1003
                    sip->next = SeqIdDup(id1);
 
1004
                } else {
 
1005
                    slp = slp1;
 
1006
                    slp1->next = slp2;
 
1007
                    sip = SeqIdDup(id1);
 
1008
                    sip->next = SeqIdDup(id2);
 
1009
                }
 
1010
                
 
1011
                sseg = StdSegNew();
 
1012
                sseg->dim = 2;
 
1013
                
 
1014
                if (sseg_head == NULL)
 
1015
                    sseg_head = sseg;
 
1016
                
 
1017
                sseg->loc = slp;
 
1018
                sseg->ids = sip;
 
1019
                
 
1020
                if (sseg_old)
 
1021
                    sseg_old->next = sseg;
 
1022
                
 
1023
                sseg_old = sseg;
 
1024
 
 
1025
                slp1 = NULL;
 
1026
                slp2 = NULL;
 
1027
            }
 
1028
 
 
1029
            first_shift = FALSE;
 
1030
 
 
1031
            /* Protein piece is empty */
 
1032
            ValNodeAddPointer(&slp1, SEQLOC_EMPTY, SeqIdDup(id1));
 
1033
            
 
1034
            /* Nucleotide scale shifted by 3, protein gapped */
 
1035
            seq_int2 = SeqIntNew();              
 
1036
            seq_int2->from = get_current_pos(&start2, curr->num*3);
 
1037
            seq_int2->to = start2 - 1;
 
1038
 
 
1039
            if(seq_int2->to >= original_length2) {
 
1040
                seq_int2->to = original_length2 -1;
 
1041
            }
 
1042
 
 
1043
            /* Transfer to DNA minus strand coordinates */
 
1044
            if(strand2 == Seq_strand_minus) {
 
1045
                int tmp_int;
 
1046
                tmp_int = seq_int2->to;
 
1047
                seq_int2->to = original_length2 - seq_int2->from - 1;
 
1048
                seq_int2->from = original_length2 - tmp_int - 1;
 
1049
            }
 
1050
 
 
1051
            seq_int2->id = SeqIdDup(id2);
 
1052
            seq_int2->strand = strand2;
 
1053
            
 
1054
            ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int2);
 
1055
            
 
1056
            seq_int1_last = NULL;
 
1057
            seq_int2_last = seq_int2; /* Will be used to adjust "to" value */
 
1058
            
 
1059
            break;
 
1060
 
 
1061
        case eGapAlignSub: /* Substitution. */
 
1062
 
 
1063
            first_shift = FALSE;
 
1064
 
 
1065
            /* Protein coordinates */
 
1066
            seq_int1 = SeqIntNew();
 
1067
            seq_int1->from =  get_current_pos(&start1, curr->num);
 
1068
            seq_int1->to = start1 - 1;
 
1069
 
 
1070
            if(seq_int1->to >= original_length1)
 
1071
                seq_int1->to = original_length1-1;
 
1072
            
 
1073
            seq_int1->id = SeqIdDup(id1);
 
1074
            seq_int1->strand = strand1;
 
1075
 
 
1076
            ValNodeAddPointer(&slp1, SEQLOC_INT, seq_int1);
 
1077
           
 
1078
            /* Nucleotide scale shifted by op_type */
 
1079
            seq_int2 = SeqIntNew();
 
1080
 
 
1081
            seq_int2->from = 
 
1082
               get_current_pos(&start2, curr->num*(Uint1)curr->op_type);
 
1083
            seq_int2->to = start2 - 1;
 
1084
 
 
1085
                /* Chop off three bases and one residue at a time.
 
1086
                        Why does this happen, seems like a bug?
 
1087
                */
 
1088
            while (seq_int2->to >= original_length2) {
 
1089
                seq_int2->to -= 3;
 
1090
                seq_int1->to--;
 
1091
            }
 
1092
 
 
1093
            /* Transfer to DNA minus strand coordinates */
 
1094
            if(strand2 == Seq_strand_minus) {
 
1095
                int tmp_int;
 
1096
                tmp_int = seq_int2->to;
 
1097
                seq_int2->to = original_length2 - seq_int2->from - 1;
 
1098
                seq_int2->from = original_length2 - tmp_int - 1;
 
1099
            }
 
1100
            
 
1101
            seq_int2->id = SeqIdDup(id2);
 
1102
            seq_int2->strand = strand2;
 
1103
 
 
1104
            ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int2);
 
1105
 
 
1106
            seq_int1_last = seq_int1; /* Will be used to adjust "to" value */
 
1107
            seq_int2_last = seq_int2; /* Will be used to adjust "to" value */
 
1108
            
 
1109
            break;
 
1110
        case eGapAlignDel2:     /* gap of two nucleotides. */
 
1111
        case eGapAlignDel1: /* Gap of one nucleotide. */
 
1112
        case eGapAlignIns1: /* Insertion of one nucleotide. */
 
1113
        case eGapAlignIns2: /* Insertion of two nucleotides. */
 
1114
 
 
1115
            if(first_shift == TRUE) { /* Second frameshift in a row */
 
1116
                /* Protein coordinates */
 
1117
                seq_int1 = SeqIntNew();
 
1118
                seq_int1->from =  get_current_pos(&start1, 1);
 
1119
                seq_int1->to = start1 - 1;
 
1120
 
 
1121
                if(seq_int1->to >= original_length1)
 
1122
                    seq_int1->to = original_length1-1;
 
1123
                
 
1124
                seq_int1->id = SeqIdDup(id1);
 
1125
                seq_int1->strand = strand1;
 
1126
                
 
1127
                ValNodeAddPointer(&slp1, SEQLOC_INT, seq_int1);
 
1128
                
 
1129
                /* Nucleotide scale shifted by op_type */
 
1130
                seq_int2 = SeqIntNew();
 
1131
 
 
1132
                seq_int2->from = 
 
1133
                   get_current_pos(&start2, (Uint1)curr->op_type);
 
1134
                seq_int2->to = start2 - 1;
 
1135
 
 
1136
                if(seq_int2->to >= original_length2) {
 
1137
                    seq_int2->to = original_length2 -1;
 
1138
                    seq_int1->to--;
 
1139
                }
 
1140
 
 
1141
                /* Transfer to DNA minus strand coordinates */
 
1142
                if(strand2 == Seq_strand_minus) {
 
1143
                    int tmp_int;
 
1144
                    tmp_int = seq_int2->to;
 
1145
                    seq_int2->to = original_length2 - seq_int2->from - 1;
 
1146
                    seq_int2->from = original_length2 - tmp_int - 1;
 
1147
                }
 
1148
            
 
1149
                seq_int2->id = SeqIdDup(id2);
 
1150
                seq_int2->strand = strand2;
 
1151
                
 
1152
                ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int2);
 
1153
 
 
1154
                seq_int1_last = seq_int1; 
 
1155
                seq_int2_last = seq_int2; 
 
1156
 
 
1157
                /* first_shift = FALSE; */
 
1158
 
 
1159
                break;
 
1160
            }
 
1161
            
 
1162
            first_shift = TRUE;
 
1163
 
 
1164
            /* If this substitution is following simple frameshift
 
1165
               we do not need to start new segment, but may continue
 
1166
               old one */
 
1167
            if(seq_int2_last != NULL) {
 
1168
                get_current_pos(&start2, curr->num*((Uint1)curr->op_type-3));
 
1169
                if(strand2 != Seq_strand_minus) {
 
1170
                    seq_int2_last->to = start2 - 1;
 
1171
                } else {
 
1172
                    /* Transfer to DNA minus strand coordinates */
 
1173
                    seq_int2_last->from = original_length2 - start2;
 
1174
                }
 
1175
 
 
1176
                /* Adjustment for multiple shifts - theoretically possible,
 
1177
                   but very unprobable */
 
1178
                if(seq_int2_last->from > seq_int2_last->to) {
 
1179
                    
 
1180
                    if(strand2 != Seq_strand_minus) {
 
1181
                        seq_int2_last->to += 3;
 
1182
                    } else {
 
1183
                        seq_int2_last->from -= 3;
 
1184
                    }
 
1185
                    
 
1186
                    if(seq_int1_last != 0)
 
1187
                        seq_int1_last++;
 
1188
                }
 
1189
 
 
1190
            } else if ((Uint1)curr->op_type > 3) {
 
1191
                /* Protein piece is empty */
 
1192
                ValNodeAddPointer(&slp1, SEQLOC_EMPTY, SeqIdDup(id1));
 
1193
                /* Simulating insertion of nucleotides */
 
1194
                seq_int2 = SeqIntNew();
 
1195
                seq_int2->from = 
 
1196
                   get_current_pos(&start2, 
 
1197
                                   curr->num*((Uint1)curr->op_type-3));
 
1198
                seq_int2->to = start2 - 1;
 
1199
                
 
1200
                if(seq_int2->to >= original_length2) {
 
1201
                    seq_int2->to = original_length2 - 1;
 
1202
                }
 
1203
 
 
1204
                /* Transfer to DNA minus strand coordinates */
 
1205
                if(strand2 == Seq_strand_minus) {
 
1206
                    int tmp_int;
 
1207
                    tmp_int = seq_int2->to;
 
1208
                    seq_int2->to = original_length2 - seq_int2->from - 1;
 
1209
                    seq_int2->from = original_length2 - tmp_int - 1;
 
1210
                }
 
1211
 
 
1212
                seq_int2->id = SeqIdDup(id2);
 
1213
                seq_int2->strand = strand2;
 
1214
                
 
1215
                ValNodeAddPointer(&slp2, SEQLOC_INT, seq_int2);
 
1216
                
 
1217
                seq_int1_last = NULL;
 
1218
                seq_int2_last = seq_int2; /* Will be used to adjust "to" value */
 
1219
                break;
 
1220
            } else {
 
1221
                continue;       /* Main loop */
 
1222
            }
 
1223
            continue;       /* Main loop */
 
1224
            /* break; */
 
1225
        default:
 
1226
            continue;       /* Main loop */
 
1227
            /* break; */
 
1228
        } 
 
1229
 
 
1230
        if (reverse) {
 
1231
            slp = slp2;
 
1232
            slp2->next = slp1;
 
1233
            sip = SeqIdDup(id2);
 
1234
            sip->next = SeqIdDup(id1);
 
1235
        } else {
 
1236
            slp = slp1;
 
1237
            slp1->next = slp2;
 
1238
            sip = SeqIdDup(id1);
 
1239
            sip->next = SeqIdDup(id2);
 
1240
        }
 
1241
 
 
1242
        sseg = StdSegNew();
 
1243
        sseg->dim = 2;
 
1244
        
 
1245
        if (sseg_head == NULL)
 
1246
            sseg_head = sseg;
 
1247
        
 
1248
        sseg->loc = slp;
 
1249
        sseg->ids = sip;
 
1250
        
 
1251
        if (sseg_old)
 
1252
            sseg_old->next = sseg;
 
1253
 
 
1254
        sseg_old = sseg;
 
1255
    }
 
1256
    sap->segs = sseg_head;
 
1257
    sap->next = NULL;
 
1258
    
 
1259
    return sap;
 
1260
}
 
1261
 
 
1262
static Int2 
 
1263
BLAST_GapInfoToSeqAlign(EBlastProgramType program_number, 
 
1264
   BlastHSPList* hsp_list, SeqIdPtr query_id, 
 
1265
   SeqIdPtr subject_id, Int4 query_length, 
 
1266
   Boolean is_ooframe, SeqAlignPtr* head_seqalign)
 
1267
{
 
1268
   Int2 status = 0;
 
1269
   BlastHSP** hsp_array;
 
1270
   SeqAlignPtr last_seqalign = NULL, seqalign = NULL;
 
1271
   Int4 index;
 
1272
 
 
1273
   *head_seqalign = NULL;
 
1274
 
 
1275
   hsp_array = hsp_list->hsp_array;
 
1276
 
 
1277
   for (index=0; index<hsp_list->hspcnt; index++) { 
 
1278
      hsp_array[index]->gap_info->original_length1 = query_length;
 
1279
      if (is_ooframe) {
 
1280
         seqalign = BLAST_OOFEditBlockToSeqAlign(program_number, 
 
1281
                       hsp_array[index]->gap_info, 
 
1282
                       query_id, subject_id);
 
1283
      } else {
 
1284
         /* The following line is needed for negative frames of translated 
 
1285
            query */
 
1286
         seqalign = GapEditBlockToSeqAlign(hsp_array[index]->gap_info, 
 
1287
                                            subject_id, query_id);
 
1288
      }
 
1289
      if (index==0) {
 
1290
         *head_seqalign = last_seqalign = seqalign;
 
1291
      } else {
 
1292
         last_seqalign->next = seqalign;
 
1293
         last_seqalign = last_seqalign->next;
 
1294
      }
 
1295
      seqalign->score = GetScoreSetFromBlastHsp(hsp_array[index]);
 
1296
   }
 
1297
 
 
1298
   return status;
 
1299
}
 
1300
 
 
1301
Int2 BLAST_ResultsToSeqAlign(EBlastProgramType program_number, 
 
1302
        BlastHSPResults* results, SeqLocPtr query_slp, 
 
1303
        ReadDBFILE* rdfp, SeqLoc* subject_slp,
 
1304
        Boolean is_gapped, Boolean is_ooframe, 
 
1305
        SeqAlignPtr* head_seqalign)
 
1306
{
 
1307
   Int4 query_index, subject_index;
 
1308
   SeqLocPtr slp = query_slp;
 
1309
   SeqIdPtr query_id, subject_id = NULL;
 
1310
   BlastHitList* hit_list;
 
1311
   BlastHSPList* hsp_list;
 
1312
   SeqAlignPtr seqalign = NULL, last_seqalign = NULL;
 
1313
   Int4 subject_length = 0;
 
1314
   SeqLoc** subject_loc_array = NULL;
 
1315
   
 
1316
   *head_seqalign = NULL;
 
1317
   if (!results)
 
1318
      return 0;
 
1319
 
 
1320
   if (!rdfp && !subject_slp)
 
1321
      return -1;
 
1322
 
 
1323
   if (!rdfp) {
 
1324
      subject_loc_array = 
 
1325
         (SeqLoc**) malloc(ValNodeLen(subject_slp)*sizeof(SeqLoc*));
 
1326
      for (slp = subject_slp, subject_index = 0; slp; slp = slp->next, ++subject_index)
 
1327
         subject_loc_array[subject_index] = slp;
 
1328
   }
 
1329
 
 
1330
   slp = query_slp;
 
1331
   for (query_index = 0; slp && query_index < results->num_queries; 
 
1332
        ++query_index, slp = slp->next) {
 
1333
      hit_list = results->hitlist_array[query_index];
 
1334
      if (!hit_list)
 
1335
         continue;
 
1336
      query_id = SeqLocId(slp);
 
1337
 
 
1338
      for (subject_index = 0; subject_index < hit_list->hsplist_count;
 
1339
           ++subject_index) {
 
1340
         hsp_list = hit_list->hsplist_array[subject_index];
 
1341
         if (!hsp_list)
 
1342
            continue;
 
1343
         if (rdfp) {
 
1344
             /* NB: The following call allocates the SeqId structure. */
 
1345
            readdb_get_descriptor(rdfp, hsp_list->oid, &subject_id, NULL);
 
1346
            subject_length = readdb_get_sequence_length(rdfp, hsp_list->oid);
 
1347
         } else {
 
1348
             /* NB: The following call does not allocate the SeqId structure,
 
1349
                but returns the existing one. */
 
1350
            subject_id = SeqLocId(subject_loc_array[hsp_list->oid]);
 
1351
            subject_length = SeqLocLen(subject_loc_array[hsp_list->oid]);
 
1352
         }
 
1353
 
 
1354
         if (is_gapped) {
 
1355
            BLAST_GapInfoToSeqAlign(program_number, hsp_list, query_id, 
 
1356
               subject_id, SeqLocLen(slp), is_ooframe, &seqalign);
 
1357
         } else {
 
1358
            BLAST_UngappedHSPToSeqAlign(program_number, hsp_list, query_id,
 
1359
               subject_id, SeqLocLen(slp), subject_length, &seqalign);
 
1360
         }                      
 
1361
 
 
1362
         /* The subject id must be deallocated only in case of a ReadDB 
 
1363
            interface */
 
1364
         if (rdfp)
 
1365
             subject_id = SeqIdSetFree(subject_id);
 
1366
 
 
1367
         if (seqalign) {
 
1368
            if (!last_seqalign) {
 
1369
               *head_seqalign = last_seqalign = seqalign;
 
1370
            } else {
 
1371
               last_seqalign->next = seqalign;
 
1372
            }
 
1373
            for ( ; last_seqalign->next; last_seqalign = last_seqalign->next);
 
1374
         }
 
1375
      }
 
1376
   }
 
1377
 
 
1378
   sfree(subject_loc_array);
 
1379
 
 
1380
   return 0;
 
1381
}
 
1382
 
 
1383
void Blast_AdjustOffsetsInSeqAlign(SeqAlign* head, SeqLoc* query_slp,
 
1384
                                   SeqLoc* subject_slp)
 
1385
{
 
1386
   SeqAlign* seqalign, *next_seqalign = NULL, *last_seqalign;
 
1387
   SeqLoc* qslp, *sslp;
 
1388
   SeqId* query_id, *subject_id;
 
1389
   
 
1390
   qslp = query_slp;
 
1391
   sslp = subject_slp;
 
1392
 
 
1393
   for (seqalign = head; seqalign; seqalign = next_seqalign) {
 
1394
      query_id = TxGetQueryIdFromSeqAlign(seqalign);
 
1395
      subject_id = TxGetSubjectIdFromSeqAlign(seqalign);
 
1396
 
 
1397
      /* Find Seq-locs corresponding to these query and subject. Start looking
 
1398
         from the previously found Seq-locs, because Seq-aligns are sorted
 
1399
         appropriately. */
 
1400
      for ( ; qslp; qslp = qslp->next) {
 
1401
         if (SeqIdComp(query_id, SeqLocId(qslp)) == SIC_YES)
 
1402
            break;
 
1403
         /* Changing query: return subject Seq-loc pointer to the beginning 
 
1404
            of the list of subject Seq-locs.*/
 
1405
         sslp = subject_slp;
 
1406
      }
 
1407
      for ( ; sslp; sslp = sslp->next) {
 
1408
         if (SeqIdComp(subject_id, SeqLocId(sslp)) == SIC_YES)
 
1409
            break;
 
1410
      }
 
1411
 
 
1412
      /* Find the first Seq-align that has either different subject or 
 
1413
         different query. */
 
1414
      last_seqalign = seqalign;
 
1415
      for (next_seqalign = seqalign->next; next_seqalign; 
 
1416
           next_seqalign = next_seqalign->next) {
 
1417
         if ((SeqIdComp(subject_id, TxGetSubjectIdFromSeqAlign(next_seqalign)) 
 
1418
              != SIC_YES) || 
 
1419
             (SeqIdComp(query_id, TxGetQueryIdFromSeqAlign(next_seqalign)) 
 
1420
              != SIC_YES)) {
 
1421
            /* Unlink the Seq-align chain */
 
1422
            last_seqalign->next = NULL;
 
1423
            break;
 
1424
         } else {
 
1425
            last_seqalign = next_seqalign;
 
1426
         }
 
1427
      }
 
1428
      AdjustOffSetsInSeqAlign(seqalign, qslp, sslp);
 
1429
      /* Reconnect the Seq-align chain */
 
1430
      last_seqalign->next = next_seqalign;
 
1431
   }
 
1432
}