← Back to branch summary

~ubuntu-branches/ubuntu/edgy/ncbi-tools6/edgy

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-19 23:28:07 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20060719232807-et3cdmcjgmnyleyx
Tags: 6.1.20060507-3ubuntu1
Re-merge with Debian

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.52 2005/04/06 23:27:53 dondosha Exp $
 
1
/* $Id: blast_seqalign.c,v 1.57 2006/03/21 15:27:58 madden Exp $
2
2
* ===========================================================================
3
3
*
4
4
*                            PUBLIC DOMAIN NOTICE
28
28
 * Conversion of BLAST results to the SeqAlign form
29
29
 */
30
30
 
31
 
static char const rcsid[] = "$Id: blast_seqalign.c,v 1.52 2005/04/06 23:27:53 dondosha Exp $";
 
31
static char const rcsid[] = "$Id: blast_seqalign.c,v 1.57 2006/03/21 15:27:58 madden Exp $";
32
32
 
33
33
#include <algo/blast/api/blast_seqalign.h>
34
34
 
35
 
extern Int4 LIBCALL HspArrayPurge PROTO((BlastHSP** hsp_array, 
36
 
                       Int4 hspcnt, Boolean clear_num));
37
35
extern SeqIdPtr GetTheSeqAlignID (SeqIdPtr seq_id);
38
36
extern ScorePtr MakeBlastScore (ScorePtr PNTR old, CharPtr scoretype, 
39
37
                                Nlm_FloatHi prob, Int4 score);
43
41
 * @{
44
42
 */
45
43
 
 
44
SBlastSeqalignArray* 
 
45
SBlastSeqalignArrayNew(Int4 size)
 
46
{
 
47
    SBlastSeqalignArray* retval = NULL;
 
48
 
 
49
    if (size <= 0)
 
50
     return retval;
 
51
 
 
52
    retval = (SBlastSeqalignArray*) malloc(sizeof(SBlastSeqalignArray));
 
53
    if (retval)
 
54
    {
 
55
        retval->num_queries = size;
 
56
        retval->array = (SeqAlign**) calloc(size, sizeof(SeqAlign*));
 
57
        if (retval->array == NULL)
 
58
        {
 
59
            sfree(retval);
 
60
            retval = NULL;
 
61
        }
 
62
    }
 
63
    return retval;
 
64
}
 
65
 
 
66
SBlastSeqalignArray* 
 
67
SBlastSeqalignArrayFree(SBlastSeqalignArray* seqalign_vec)
 
68
{
 
69
   Int4 index;
 
70
 
 
71
   if (seqalign_vec == NULL)
 
72
      return NULL;
 
73
 
 
74
   for (index=0; index<seqalign_vec->num_queries; index++)
 
75
   {
 
76
        SeqAlignSetFree(seqalign_vec->array[index]); 
 
77
   }
 
78
   sfree(seqalign_vec->array);
 
79
   sfree(seqalign_vec);
 
80
   return NULL;
 
81
}
 
82
 
46
83
/** Creates a score set corresponding to one HSP.
47
84
 * @param hsp HSP structure [in]
48
85
 * @return Score set for this HSP.
83
120
   
84
121
   if (hsp->num_ident > 0)
85
122
      MakeBlastScore(&score_set, "num_ident", 0.0, hsp->num_ident);
86
 
   
 
123
 
 
124
   if (hsp->comp_adjustment_method > 0) { 
 
125
       MakeBlastScore(&score_set, "comp_adjustment_method",0.0,
 
126
                      hsp->comp_adjustment_method);
 
127
   }
87
128
   return score_set;
88
129
}
89
130
 
383
424
}
384
425
 
385
426
Int2 
386
 
GapCollectDataForSeqalign(BlastHSP* hsp, GapEditScript* curr_in, Int4 numseg,
 
427
GapCollectDataForSeqalign(BlastHSP* hsp, GapEditScript* esp, Int4 first, Int4 number,
387
428
                          Int4 query_length, Int4 subject_length,
388
429
                          Boolean translate1, Boolean translate2,
389
430
                          Int4** start_out, Int4** length_out, 
390
431
                          Uint1** strands_out, Int4* start1, Int4* start2)
391
432
{
392
 
    GapEditScript* curr;
393
433
    Int2 frame1, frame2;
394
434
    Int4 begin1, begin2, index, length1, length2;
395
435
    Int4 original_length1, original_length2, i;
423
463
    else
424
464
        strand2 = Seq_strand_unknown; 
425
465
 
426
 
    start = (Int4 *) calloc((2*numseg+1), sizeof(Int4));
427
 
    length = (Int4 *) calloc((numseg+1), sizeof(Int4));
428
 
    strands = (Uint1 *) calloc((2*numseg+1), sizeof(Uint1));
 
466
    start = (Int4 *) calloc((2*number+1), sizeof(Int4));
 
467
    length = (Int4 *) calloc((number+1), sizeof(Int4));
 
468
    strands = (Uint1 *) calloc((2*number+1), sizeof(Uint1));
429
469
 
430
470
    index=0;
431
 
    for (i = 0, curr=curr_in; curr && i < numseg; curr=curr->next, i++) {
432
 
        switch(curr->op_type) {
 
471
    for (i=first; i<number; i++)
 
472
    {
 
473
        switch(esp->op_type[i]) {
433
474
        case eGapAlignDecline:
434
475
        case eGapAlignSub:
435
476
            if (strand1 != Seq_strand_minus) {
436
477
                if(translate1 == FALSE)
437
 
                    begin1 = s_GetCurrentPos(start1, curr->num);
 
478
                    begin1 = s_GetCurrentPos(start1, esp->num[i]);
438
479
                else
439
 
                    begin1 = frame1 - 1 + CODON_LENGTH*s_GetCurrentPos(start1, curr->num);
 
480
                    begin1 = frame1 - 1 + CODON_LENGTH*s_GetCurrentPos(start1, esp->num[i]);
440
481
            } else {
441
482
                if(translate1 == FALSE)
442
 
                    begin1 = length1 - s_GetCurrentPos(start1, curr->num) - curr->num;
 
483
                    begin1 = length1 - s_GetCurrentPos(start1, esp->num[i]) - esp->num[i];
443
484
                else
444
 
                    begin1 = original_length1 - CODON_LENGTH*(s_GetCurrentPos(start1, curr->num)+curr->num) + frame1 + 1;
 
485
                    begin1 = original_length1 - CODON_LENGTH*(s_GetCurrentPos(start1, esp->num[i])+esp->num[i]) + frame1 + 1;
445
486
            }
446
487
                                        
447
488
            if (strand2 != Seq_strand_minus) {
448
489
                if(translate2 == FALSE)
449
 
                    begin2 = s_GetCurrentPos(start2, curr->num);
 
490
                    begin2 = s_GetCurrentPos(start2, esp->num[i]);
450
491
                else
451
 
                    begin2 = frame2 - 1 + CODON_LENGTH*s_GetCurrentPos(start2, curr->num);
 
492
                    begin2 = frame2 - 1 + CODON_LENGTH*s_GetCurrentPos(start2, esp->num[i]);
452
493
            } else {
453
494
                if(translate2 == FALSE)
454
 
                    begin2 = length2 - s_GetCurrentPos(start2, curr->num) - curr->num;
 
495
                    begin2 = length2 - s_GetCurrentPos(start2, esp->num[i]) - esp->num[i];
455
496
                else
456
 
                    begin2 = original_length2 - CODON_LENGTH*(s_GetCurrentPos(start2, curr->num)+curr->num) + frame2 + 1;
 
497
                    begin2 = original_length2 - CODON_LENGTH*(s_GetCurrentPos(start2, esp->num[i])+esp->num[i]) + frame2 + 1;
457
498
            }
458
499
            
459
500
            strands[2*index] = strand1;
467
508
            begin1 = -1;
468
509
            if (strand2 != Seq_strand_minus) {
469
510
                if(translate2 == FALSE)
470
 
                    begin2 = s_GetCurrentPos(start2, curr->num);
 
511
                    begin2 = s_GetCurrentPos(start2, esp->num[i]);
471
512
                else
472
 
                    begin2 = frame2 - 1 + CODON_LENGTH*s_GetCurrentPos(start2, curr->num);
 
513
                    begin2 = frame2 - 1 + CODON_LENGTH*s_GetCurrentPos(start2, esp->num[i]);
473
514
            } else {
474
515
                if(translate2 == FALSE)
475
 
                    begin2 = length2 - s_GetCurrentPos(start2, curr->num) - curr->num;
 
516
                    begin2 = length2 - s_GetCurrentPos(start2, esp->num[i]) - esp->num[i];
476
517
                else
477
 
                    begin2 = original_length2 - CODON_LENGTH*(s_GetCurrentPos(start2, curr->num)+curr->num) + frame2 + 1;
 
518
                    begin2 = original_length2 - CODON_LENGTH*(s_GetCurrentPos(start2, esp->num[i])+esp->num[i]) + frame2 + 1;
478
519
            }
479
520
            
480
521
            if (index > 0)
490
531
        case eGapAlignIns:
491
532
            if (strand1 != Seq_strand_minus) {
492
533
                if(translate1 == FALSE)
493
 
                    begin1 = s_GetCurrentPos(start1, curr->num);
 
534
                    begin1 = s_GetCurrentPos(start1, esp->num[i]);
494
535
                else
495
 
                    begin1 = frame1 - 1 + CODON_LENGTH*s_GetCurrentPos(start1, curr->num);
 
536
                    begin1 = frame1 - 1 + CODON_LENGTH*s_GetCurrentPos(start1, esp->num[i]);
496
537
            } else {
497
538
                if(translate1 == FALSE)
498
 
                    begin1 = length1 - s_GetCurrentPos(start1, curr->num) - curr->num;
 
539
                    begin1 = length1 - s_GetCurrentPos(start1, esp->num[i]) - esp->num[i];
499
540
                else
500
 
                    begin1 = original_length1 - CODON_LENGTH*(s_GetCurrentPos(start1, curr->num)+curr->num) + frame1 + 1;
 
541
                    begin1 = original_length1 - CODON_LENGTH*(s_GetCurrentPos(start1, esp->num[i])+esp->num[i]) + frame1 + 1;
501
542
            }
502
543
            begin2 = -1;
503
544
            strands[2*index] = strand1;
512
553
        default:
513
554
            break;
514
555
        }
515
 
        length[index] = curr->num;
 
556
        length[index] = esp->num[i];
516
557
        index++;
517
558
    }    
518
559
 
538
579
static void 
539
580
s_GapCorrectUASequence(BlastHSP* hsp)
540
581
{
541
 
    GapEditScript* curr,* curr_last,* curr_last2;
542
 
    Boolean last_indel;
543
 
 
544
 
    last_indel = FALSE;
545
 
    curr_last = NULL;
546
 
    curr_last2 = NULL;
547
 
 
548
 
    for (curr=hsp->gap_info; curr; curr = curr->next) {
549
 
 
550
 
        if(curr->op_type == eGapAlignDecline && last_indel == TRUE) {
 
582
    GapEditScript* esp = hsp->gap_info;
 
583
    int index;
 
584
 
 
585
    for (index=0; index<esp->size; index++)
 
586
    {
 
587
        // if GAPALIGN_DECLINE immediately follows an insertion or deletion
 
588
        if (index > 0 && esp->op_type[index] == eGapAlignDecline &&
 
589
           (esp->op_type[index-1] == eGapAlignIns || esp->op_type[index-1] == eGapAlignDel))
 
590
        {
551
591
            /* This is invalid condition and regions should be
552
592
               exchanged */
553
 
 
554
 
            if(curr_last2 != NULL)
555
 
                curr_last2->next = curr;
556
 
            else
557
 
                hsp->gap_info = curr; /* First element in a list */
558
 
            
559
 
            curr_last->next = curr->next;
560
 
            curr->next = curr_last;
561
 
        }
562
 
        
563
 
        last_indel = FALSE;
564
 
        
565
 
        if(curr->op_type == eGapAlignIns || curr->op_type == eGapAlignDel) {
566
 
            last_indel = TRUE;
567
 
            curr_last2 = curr_last;
568
 
        }
569
 
 
570
 
        curr_last = curr;
 
593
            int temp_num = esp->num[index];
 
594
            EGapAlignOpType temp_op = esp->op_type[index];
 
595
 
 
596
            esp->num[index] = esp->num[index-1];
 
597
            esp->op_type[index] = esp->op_type[index-1];
 
598
            esp->num[index-1] = temp_num;
 
599
            esp->op_type[index-1] = temp_op;
 
600
        }
571
601
    }
572
602
    return;
573
603
}
699
729
                     Int4 query_length, Int4 subject_length)
700
730
 
701
731
{
702
 
    GapEditScript* curr,* curr_last;
703
 
    Int4 numseg, start1, start2;
 
732
    GapEditScript* esp;
 
733
    Int4 start1, start2;
704
734
    Int4* length,* start;
705
735
    Uint1* strands;
706
736
    Boolean is_disc_align = FALSE;
707
737
    SeqAlignPtr sap, sap_disc, sap_head, sap_tail;
708
 
    Boolean skip_region = FALSE;
709
738
    Boolean translate1, translate2;
710
 
 
711
 
    is_disc_align = FALSE; numseg = 0;
712
 
 
713
 
    for (curr=hsp->gap_info; curr; curr = curr->next) {
714
 
        numseg++;
715
 
        if(/*edit_block->discontinuous && */curr->op_type == eGapAlignDecline)
 
739
    int index;
 
740
 
 
741
    is_disc_align = FALSE;
 
742
 
 
743
    esp = hsp->gap_info;
 
744
    for (index=0; index<esp->size; index++)
 
745
    {
 
746
        if(esp->op_type[index] == eGapAlignDecline)
 
747
        {
716
748
           is_disc_align = TRUE;
 
749
           break;
 
750
        }
717
751
    }
718
752
    
719
753
    start1 = hsp->query.offset;
730
764
           strand, translate, reverse etc. Real data is taken starting
731
765
           from "curr" and taken only "numseg" segments */
732
766
        
733
 
        GapCollectDataForSeqalign(hsp, hsp->gap_info, numseg, query_length,
 
767
        GapCollectDataForSeqalign(hsp, hsp->gap_info, 0, esp->size, query_length,
734
768
                                  subject_length, translate1, translate2,
735
769
                                  &start, &length, &strands, &start1, &start2);
736
770
        
737
771
        /* Result of this function will be either den-seg or Std-seg
738
772
           depending on translation options */
739
773
        sap = s_GapMakeSeqAlign(query_id, subject_id, translate1, translate2, 
740
 
                                numseg, length, start, strands);
 
774
                                esp->size, length, start, strands);
741
775
    } else {
742
776
 
743
777
        /* By request of Steven Altschul - we need to have 
752
786
        sap_disc->type = SAT_PARTIAL; /* ordered segments, over part of seq */
753
787
        sap_disc->segtype = SAS_DISC; /* discontinuous alignment */
754
788
        
755
 
        curr_last = hsp->gap_info; 
756
789
        sap_head = NULL; sap_tail = NULL;
757
 
        while(curr_last != NULL) {
758
 
            skip_region = FALSE;                        
759
 
            for (numseg = 0, curr = curr_last; 
760
 
                 curr; curr = curr->next, numseg++) {
 
790
        for (index=0; index<esp->size; index++)
 
791
        {
 
792
            int numseg=0;
 
793
            Boolean skip_region = FALSE;
 
794
            int index2 = index;
 
795
            int first = index;
 
796
            for (index2=first; index2<esp->size; index2++, numseg++) {
761
797
 
762
 
                if(curr->op_type == eGapAlignDecline) {
 
798
                if(esp->op_type[index2] == eGapAlignDecline) {
763
799
                    if(numseg != 0) { /* End of aligned area */
764
800
                        break;
765
801
                    } else {
766
 
                        while(curr && curr->op_type == eGapAlignDecline) {
 
802
                        while (index2<esp->size && esp->op_type[index2] == eGapAlignDecline) {
767
803
                            numseg++;
768
 
                            curr = curr->next; 
 
804
                            index2++;
769
805
                        }
770
806
                        skip_region = TRUE;                        
771
807
                        break;
772
808
                    }
773
809
                }
774
810
            }
775
 
 
776
 
            GapCollectDataForSeqalign(hsp, curr_last, numseg, query_length,
 
811
            
 
812
 
 
813
            if(!skip_region) {            
 
814
 
 
815
               GapCollectDataForSeqalign(hsp, esp, first, numseg, query_length,
777
816
                                      subject_length, translate1, translate2,
778
817
                                      &start, &length, &strands, &start1, 
779
818
                                      &start2);
780
819
            
781
 
            if(!skip_region) {            
782
820
                sap = 
783
821
                    s_GapMakeSeqAlign(query_id, subject_id, translate1, 
784
822
                                      translate2, numseg, length, start, 
792
830
                    sap_tail = sap;
793
831
                }
794
832
            }
795
 
            
796
 
            curr_last = curr;
797
833
        }
798
834
        sap_disc->segs = sap_head;
799
835
        sap = sap_disc;
818
854
                         Int4 query_length, Int4 subject_length)
819
855
{
820
856
    Boolean reverse = FALSE;
821
 
    GapEditScript* curr,* esp;
 
857
    GapEditScript* esp;
822
858
    Int2 frame1, frame2;
823
859
    Int4 start1, start2;
824
860
    Int4 original_length1, original_length2;
830
866
    StdSeg* sseg,* sseg_head,* sseg_old;
831
867
    Uint1 strand1, strand2;
832
868
    Boolean first_shift;
 
869
    int index;
833
870
 
834
871
    if (program == eBlastTypeBlastx) {
835
872
       reverse = TRUE;
866
903
    else
867
904
        strand2 = Seq_strand_unknown;
868
905
    
869
 
    esp = hsp->gap_info;
870
906
    
871
907
    sap = SeqAlignNew();
872
908
    
879
915
 
880
916
    first_shift = FALSE;
881
917
 
882
 
    for (curr=esp; curr; curr=curr->next) {
 
918
    esp = hsp->gap_info;
883
919
 
 
920
    for (index=0; index<esp->size; index++)
 
921
    {
884
922
        slp1 = NULL;
885
923
        slp2 = NULL;
886
924
        
887
 
        switch (curr->op_type) {
 
925
        switch (esp->op_type[index]) {
888
926
        case eGapAlignDel: /* deletion of three nucleotides. */
889
927
            
890
928
            first_shift = FALSE;
891
929
 
892
930
            seq_int1 = SeqIntNew();
893
 
            seq_int1->from = s_GetCurrentPos(&start1, curr->num);
 
931
            seq_int1->from = s_GetCurrentPos(&start1, esp->num[index]);
894
932
            seq_int1->to = start1 - 1;            
895
933
 
896
934
            if(seq_int1->to >= original_length1)
995
1033
            
996
1034
            /* Nucleotide scale shifted by 3, protein gapped */
997
1035
            seq_int2 = SeqIntNew();              
998
 
            seq_int2->from = s_GetCurrentPos(&start2, curr->num*3);
 
1036
            seq_int2->from = s_GetCurrentPos(&start2, esp->num[index]*3);
999
1037
            seq_int2->to = start2 - 1;
1000
1038
 
1001
1039
            if(seq_int2->to >= original_length2) {
1026
1064
 
1027
1065
            /* Protein coordinates */
1028
1066
            seq_int1 = SeqIntNew();
1029
 
            seq_int1->from =  s_GetCurrentPos(&start1, curr->num);
 
1067
            seq_int1->from =  s_GetCurrentPos(&start1, esp->num[index]);
1030
1068
            seq_int1->to = start1 - 1;
1031
1069
 
1032
1070
            if(seq_int1->to >= original_length1)
1041
1079
            seq_int2 = SeqIntNew();
1042
1080
 
1043
1081
            seq_int2->from = 
1044
 
               s_GetCurrentPos(&start2, curr->num*(Uint1)curr->op_type);
 
1082
               s_GetCurrentPos(&start2, esp->num[index]*(Uint1)esp->op_type[index]);
1045
1083
            seq_int2->to = start2 - 1;
1046
1084
 
1047
1085
                /* Chop off three bases and one residue at a time.
1092
1130
                seq_int2 = SeqIntNew();
1093
1131
 
1094
1132
                seq_int2->from = 
1095
 
                   s_GetCurrentPos(&start2, (Uint1)curr->op_type);
 
1133
                   s_GetCurrentPos(&start2, (Uint1)esp->op_type[index]);
1096
1134
                seq_int2->to = start2 - 1;
1097
1135
 
1098
1136
                if(seq_int2->to >= original_length2) {
1127
1165
               we do not need to start new segment, but may continue
1128
1166
               old one */
1129
1167
            if(seq_int2_last != NULL) {
1130
 
                s_GetCurrentPos(&start2, curr->num*((Uint1)curr->op_type-3));
 
1168
                s_GetCurrentPos(&start2, esp->num[index]*((Uint1)esp->op_type[index]-3));
1131
1169
                if(strand2 != Seq_strand_minus) {
1132
1170
                    seq_int2_last->to = start2 - 1;
1133
1171
                } else {
1149
1187
                        seq_int1_last++;
1150
1188
                }
1151
1189
 
1152
 
            } else if ((Uint1)curr->op_type > 3) {
 
1190
            } else if ((Uint1)esp->op_type[index] > 3) {
1153
1191
                /* Protein piece is empty */
1154
1192
                ValNodeAddPointer(&slp1, SEQLOC_EMPTY, SeqIdDup(id1));
1155
1193
                /* Simulating insertion of nucleotides */
1156
1194
                seq_int2 = SeqIntNew();
1157
1195
                seq_int2->from = 
1158
1196
                   s_GetCurrentPos(&start2, 
1159
 
                                   curr->num*((Uint1)curr->op_type-3));
 
1197
                                   esp->num[index]*((Uint1)esp->op_type[index]-3));
1160
1198
                seq_int2->to = start2 - 1;
1161
1199
                
1162
1200
                if(seq_int2->to >= original_length2) {
1273
1311
   return status;
1274
1312
}
1275
1313
 
1276
 
/** Adjusts the offsets in a Seq-align list produced by a BLAST search
1277
 
 * of multiple queries against multiple subjects. Seq-aligns in a list are
1278
 
 * assumed to be sorted by query, and then by subject, i.e. all Seq-aligns for
1279
 
 * the first query are at the beginning of the list, then all for the second
1280
 
 * query, etc. Within a list for any single query, all Seq-aligns for one 
1281
 
 * subject are grouped together. The order of queries in the Seq-align list 
1282
 
 * must be the same as in the query Seq-loc list; the order of subjects
1283
 
 * within a Seq-align list for a given query is the same as in the subject
1284
 
 * Seq-loc list. Some or all queries or subjects might be skipped in the 
1285
 
 * Seq-align list.
1286
 
 * @param head List of Seq-aligns [in|out]
1287
 
 * @param query_slp List of query locations [in]
1288
 
 * @param subject_slp List of subject locations [in]
1289
 
 */
1290
 
static void 
1291
 
s_AdjustOffsetsInSeqAlign(SeqAlign* head, SeqLoc* query_slp,
1292
 
                          SeqLoc* subject_slp)
1293
 
{
1294
 
   SeqAlign* seqalign, *next_seqalign = NULL, *last_seqalign;
1295
 
   SeqLoc* qslp, *sslp;
1296
 
   SeqId* query_id, *subject_id;
1297
 
   
1298
 
   qslp = query_slp;
1299
 
   sslp = subject_slp;
1300
 
 
1301
 
   for (seqalign = head; seqalign; seqalign = next_seqalign) {
1302
 
      query_id = TxGetQueryIdFromSeqAlign(seqalign);
1303
 
      subject_id = TxGetSubjectIdFromSeqAlign(seqalign);
1304
 
 
1305
 
      /* Find Seq-locs corresponding to these query and subject. Start looking
1306
 
         from the previously found Seq-locs, because Seq-aligns are sorted
1307
 
         appropriately. */
1308
 
      for ( ; qslp; qslp = qslp->next) {
1309
 
         if (SeqIdComp(query_id, SeqLocId(qslp)) == SIC_YES)
1310
 
            break;
1311
 
         /* Changing query: return subject Seq-loc pointer to the beginning 
1312
 
            of the list of subject Seq-locs.*/
1313
 
         sslp = subject_slp;
1314
 
      }
1315
 
      for ( ; sslp; sslp = sslp->next) {
1316
 
         if (SeqIdComp(subject_id, SeqLocId(sslp)) == SIC_YES)
1317
 
            break;
1318
 
      }
1319
 
 
1320
 
      /* Find the first Seq-align that has either different subject or 
1321
 
         different query. */
1322
 
      last_seqalign = seqalign;
1323
 
      for (next_seqalign = seqalign->next; next_seqalign; 
1324
 
           next_seqalign = next_seqalign->next) {
1325
 
         if ((SeqIdComp(subject_id, TxGetSubjectIdFromSeqAlign(next_seqalign)) 
1326
 
              != SIC_YES) || 
1327
 
             (SeqIdComp(query_id, TxGetQueryIdFromSeqAlign(next_seqalign)) 
1328
 
              != SIC_YES)) {
1329
 
            /* Unlink the Seq-align chain */
1330
 
            last_seqalign->next = NULL;
1331
 
            break;
1332
 
         } else {
1333
 
            last_seqalign = next_seqalign;
1334
 
         }
1335
 
      }
1336
 
      AdjustOffSetsInSeqAlign(seqalign, qslp, sslp);
1337
 
      /* Reconnect the Seq-align chain */
1338
 
      last_seqalign->next = next_seqalign;
1339
 
   }
1340
 
}
1341
 
 
1342
1314
Int2 BLAST_ResultsToSeqAlign(EBlastProgramType program_number, 
1343
 
        BlastHSPResults* results, SeqLocPtr query_slp, 
 
1315
        BlastHSPResults** results_ptr, SeqLocPtr query_slp, 
1344
1316
        ReadDBFILE* rdfp, SeqLoc* subject_slp,
1345
1317
        Boolean is_gapped, Boolean is_ooframe, 
1346
 
        SeqAlignPtr* head_seqalign)
 
1318
        SBlastSeqalignArray* *seqalign_arr)
1347
1319
{
1348
1320
   Int4 query_index, subject_index;
1349
1321
   SeqLocPtr slp = query_slp;
1350
1322
   SeqIdPtr query_id, subject_id = NULL;
1351
 
   BlastHitList* hit_list;
1352
 
   BlastHSPList* hsp_list;
1353
 
   SeqAlignPtr seqalign = NULL, last_seqalign = NULL;
1354
 
   Int4 subject_length = 0;
1355
1323
   SeqLoc** subject_loc_array = NULL;
1356
 
   
1357
 
   *head_seqalign = NULL;
1358
 
   if (!results)
 
1324
   BlastHSPResults* results = NULL;
 
1325
   
 
1326
   if (results_ptr == NULL)
 
1327
      return 0;
 
1328
 
 
1329
   results = *results_ptr;
 
1330
   
 
1331
   if (!results || results->num_queries <= 0)
1359
1332
      return 0;
1360
1333
 
1361
1334
   if (!rdfp && !subject_slp)
1362
1335
      return -1;
1363
1336
 
 
1337
   *seqalign_arr = SBlastSeqalignArrayNew(results->num_queries);
 
1338
   if (*seqalign_arr == NULL)
 
1339
      return -1;
 
1340
   
 
1341
 
1364
1342
   if (!rdfp) {
1365
1343
      subject_loc_array = 
1366
1344
         (SeqLoc**) malloc(ValNodeLen(subject_slp)*sizeof(SeqLoc*));
1371
1349
   slp = query_slp;
1372
1350
   for (query_index = 0; slp && query_index < results->num_queries; 
1373
1351
        ++query_index, slp = slp->next) {
1374
 
      hit_list = results->hitlist_array[query_index];
 
1352
      SeqAlignPtr head_seqalign = NULL, last_seqalign = NULL;
 
1353
      BlastHitList* hit_list = results->hitlist_array[query_index];
1375
1354
      if (!hit_list)
1376
1355
         continue;
 
1356
 
1377
1357
      query_id = SeqLocId(slp);
1378
1358
 
1379
1359
      for (subject_index = 0; subject_index < hit_list->hsplist_count;
1380
1360
           ++subject_index) {
1381
 
         hsp_list = hit_list->hsplist_array[subject_index];
 
1361
         SeqAlignPtr seqalign = NULL;
 
1362
         Int4 subject_length = 0;
 
1363
         BlastHSPList* hsp_list = hit_list->hsplist_array[subject_index];
1382
1364
         if (!hsp_list)
1383
1365
            continue;
1384
1366
 
1385
 
         // Sort HSPs with e-values as first priority and scores as 
1386
 
         // tie-breakers, since that is the order we want to see them in 
1387
 
         // in Seq-aligns.
 
1367
         /* Sort HSPs with e-values as first priority and scores as 
 
1368
            tie-breakers, since that is the order we want to see them in 
 
1369
            in Seq-aligns. */
1388
1370
         Blast_HSPListSortByEvalue(hsp_list);
1389
1371
 
1390
1372
         if (rdfp) {
1408
1390
                                        subject_length, &seqalign);
1409
1391
         }                      
1410
1392
 
 
1393
         if (seqalign)
 
1394
         {
 
1395
             SeqLocPtr subject_loc = NULL;
 
1396
             if (subject_loc_array)
 
1397
                subject_loc = subject_loc_array[hsp_list->oid];
 
1398
             AdjustOffSetsInSeqAlign(seqalign, slp, subject_loc); 
 
1399
         }
 
1400
 
1411
1401
         /* The subject id must be deallocated only in case of a ReadDB 
1412
1402
            interface */
1413
1403
         if (rdfp)
1415
1405
 
1416
1406
         if (seqalign) {
1417
1407
            if (!last_seqalign) {
1418
 
               *head_seqalign = last_seqalign = seqalign;
 
1408
               head_seqalign = last_seqalign = seqalign;
1419
1409
            } else {
1420
1410
               last_seqalign->next = seqalign;
1421
1411
            }
1422
1412
            for ( ; last_seqalign->next; last_seqalign = last_seqalign->next);
1423
1413
         }
1424
1414
      }
 
1415
      (*seqalign_arr)->array[query_index] = head_seqalign;
 
1416
      results->hitlist_array[query_index] = Blast_HitListFree(results->hitlist_array[query_index]);
1425
1417
   }
1426
1418
 
 
1419
   results = Blast_HSPResultsFree(results);
 
1420
   *results_ptr = NULL;
1427
1421
   sfree(subject_loc_array);
1428
1422
 
1429
 
   s_AdjustOffsetsInSeqAlign(*head_seqalign, query_slp, subject_slp);
1430
 
 
1431
1423
   return 0;
1432
1424
}
1433
1425
/* @} */