← Back to branch summary

~ubuntu-branches/ubuntu/warty/ncbi-tools6/warty

« back to all changes in this revision

Viewing changes to tools/salign.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2002-04-04 22:13:09 UTC
  • Revision ID: james.westby@ubuntu.com-20020404221309-vfze028rfnlrldct
Tags: upstream-6.1.20011220a
Import upstream version 6.1.20011220a

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/salign.c
 
1
/*   salign.c
 
2
* ===========================================================================
 
3
*
 
4
*                            PUBLIC DOMAIN NOTICE
 
5
*            National Center for Biotechnology Information (NCBI)
 
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 official 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 do not place any restriction on its use or reproduction.
 
13
*  We would, however, appreciate having the NCBI and the author cited in
 
14
*  any work or product based on this material
 
15
*
 
16
*  Although all reasonable efforts have been taken to ensure the accuracy
 
17
*  and reliability of the software and data, the NLM and the U.S.
 
18
*  Government do not and cannot warrant the performance or results that
 
19
*  may be obtained by using this software or data. The NLM and the U.S.
 
20
*  Government disclaim all warranties, express or implied, including
 
21
*  warranties of performance, merchantability or fitness for any particular
 
22
*  purpose.
 
23
*
 
24
* ===========================================================================
 
25
*
 
26
* File Name:  salign.c
 
27
*
 
28
* Author:  Colombe Chappey 
 
29
*
 
30
* Version Creation Date:   8/22/99
 
31
*
 
32
* File Description: 
 
33
*
 
34
* Modifications:  
 
35
* --------------------------------------------------------------------------
 
36
* Date     Name        Description of modification
 
37
* -------  ----------  -----------------------------------------------------
 
38
*
 
39
*
 
40
* ==========================================================================
 
41
*/
 
42
 
 
43
#include <salign.h>
 
44
#include <ncbi.h>
 
45
#include <objalign.h>
 
46
#include <bandalgn.h>
 
47
#include <blast.h>
 
48
#include <salutil.h>
 
49
#include <salsap.h>
 
50
#include <salstruc.h>
 
51
#include <salpedit.h>
 
52
 
 
53
#ifdef SALSA_DEBUG
 
54
#include <simutil.h>
 
55
#include <toasn3.h>
 
56
#include <utilpub.h>
 
57
#include <tfuns.h>
 
58
#endif
 
59
 
 
60
 
 
61
#define BAND_LIMIT 100
 
62
 
 
63
typedef struct nodehit
 
64
{
 
65
  Int2        index;
 
66
  Boolean     open;
 
67
  SeqAlignPtr salp;
 
68
  Int4        score;
 
69
  Nlm_FloatHi bit_score;
 
70
  Nlm_FloatHi evalue;
 
71
  struct nodehit PNTR     child;
 
72
  struct nodehit PNTR     next;
 
73
} NodeHit, PNTR NodeHitPtr;
 
74
 
 
75
static ValNodePtr traverseGraph (NodeHitPtr head, NodeHitPtr node, Int2 total, Int2 level, Int2Ptr tab, ValNodePtr vnp, Boolean * change);
 
76
 
 
77
NLM_EXTERN MashPtr MashNew (Boolean is_prot)
 
78
{
 
79
  MashPtr msp;
 
80
  
 
81
    if((msp = (MashPtr)MemNew(sizeof(Mash)))==NULL) {
 
82
      ErrPostEx(SEV_WARNING, 0, 0, "MemNew returned NULL");
 
83
      return NULL;
 
84
    }
 
85
     msp->is_prot = is_prot;
 
86
     msp->band_method = 2;
 
87
     msp->multidim = FALSE;
 
88
     if (msp->is_prot) 
 
89
     {
 
90
        msp->reward = 1;
 
91
        msp->penalty = -2;
 
92
        msp->gap_open = 11;
 
93
        msp->gap_extend = 1;
 
94
        msp->matrixname="BLOSUM62";
 
95
        msp->wordsize = 3;
 
96
        msp->gap_x_dropoff = 15;
 
97
        msp->gap_x_dropoff_final = 500;   /****25****/
 
98
        msp->filter = FILTER_SEG;
 
99
     } 
 
100
     else {
 
101
        msp->reward = 2;
 
102
        msp->penalty = -3;
 
103
        msp->gap_open = 10;
 
104
        msp->gap_extend = 2;
 
105
        msp->matrixname=NULL;
 
106
        msp->wordsize = 11;
 
107
        msp->gap_x_dropoff = 50;
 
108
        msp->gap_x_dropoff_final = 50;
 
109
        msp->filter = FILTER_DUST;
 
110
     }
 
111
     msp->translate_prot = FALSE;
 
112
     msp->transalp = NULL;
 
113
     msp->use_gapped_blast = TRUE;
 
114
     msp->lg1_ext = 0; 
 
115
     msp->rg1_ext =0; 
 
116
     msp->lg2_ext =0; 
 
117
     msp->rg2_ext =0; 
 
118
     msp->lg1_open =0; 
 
119
     msp->lg2_open =0; 
 
120
     msp->rg1_open =0; 
 
121
     msp->rg2_open =0; 
 
122
     msp->blast_threshold = 50;
 
123
     msp->choice_blastfilter = 2;      /* 1 */
 
124
     msp->splicing = TRUE;              /*FALSE;*/
 
125
     msp->map_align = FALSE;
 
126
     msp->align_ends = TRUE;
 
127
     return msp;
 
128
}
 
129
 
 
130
static Int4Ptr PNTR LIBCALL BlastStyleMatDestruct(Int4Ptr PNTR matrix){
 
131
   MemFree(matrix);
 
132
   matrix=NULL;
 
133
   return matrix;
 
134
}
 
135
 
 
136
static GlobalBandStructPtr LIBCALL DestructBandStruct(GlobalBandStructPtr gbsp){
 
137
   gbsp->matrix = BlastStyleMatDestruct(gbsp->matrix);
 
138
   if(gbsp->seq2!=NULL) MemFree(gbsp->seq2);
 
139
   if(gbsp->seq1!=NULL) MemFree(gbsp->seq1);
 
140
   gbsp->seq1=NULL;
 
141
   gbsp->seq2=NULL;
 
142
   gbsp=GlobalBandStructDelete(gbsp);
 
143
   return gbsp;
 
144
}
 
145
 
 
146
/* BANDALIGN */
 
147
static GlobalBandStructPtr CC_CreatBandStruct(SeqLocPtr slp1, SeqLocPtr slp2, MashPtr msp)
 
148
{
 
149
   GlobalBandStructPtr gbsp;
 
150
   Boolean is_prot;
 
151
 
 
152
   is_prot = (Boolean)(msp->is_prot || msp->translate_prot);
 
153
   gbsp = CreatBandStruct(slp1, slp2, NULL,(Boolean) is_prot, (Int2) msp->band_method);
 
154
   if (( ChangeGlobalBandMatrix(gbsp, is_prot, msp->matrixname,(Int4) msp->penalty, (Int4)msp->reward)) != 0) {
 
155
      gbsp = GlobalBandStructDelete (gbsp);
 
156
      return NULL;
 
157
   }
 
158
   SetGlobaltOptions(gbsp,msp->lg1_ext,msp->rg1_ext,
 
159
                          msp->lg2_ext, msp->rg2_ext, 
 
160
                          msp->lg1_open, msp->lg2_open, 
 
161
                          msp->rg1_open, msp->rg2_open,
 
162
                          (Int2)msp->gap_open, (Int2) msp->gap_extend);
 
163
   return gbsp;
 
164
}
 
165
/********************************************************/
 
166
static SeqAlignPtr BandAlignTwoSeqLocsFunc (SeqLocPtr slp1, SeqLocPtr slp2, MashPtr msp)
 
167
{
 
168
  GlobalBandStructPtr gbsp;
 
169
  SeqAlignPtr         newalign = NULL;
 
170
  SeqIntPtr sit;
 
171
  Boolean   is_prot;
 
172
 
 
173
  if (slp1==NULL || slp2==NULL)
 
174
     return NULL;
 
175
  gbsp = CC_CreatBandStruct(slp1, slp2, msp);
 
176
  if (gbsp!=NULL) 
 
177
  {
 
178
    is_prot = (Boolean)(msp->is_prot || msp->translate_prot);
 
179
    gbsp->options->low = -(Int4) SeqLocLen(slp1);
 
180
    gbsp->options->up =(Int4) SeqLocLen(slp2);
 
181
    gbsp->options->start_diag = gbsp->options->low;
 
182
    gbsp->options->width = gbsp->options->up-gbsp->options->low + 1;
 
183
/*     
 
184
    if (gbsp->options->width>2*BAND_LIMIT) {
 
185
      SetLowUpFromBlast(gbsp->options, is_prot, (Int2)0, (Int2)30,slp1, slp2);
 
186
    }
 
187
    else 
 
188
    {
 
189
       gbsp->search_type = (Uint1) G_BAND_QGAP; 
 
190
    }
 
191
*/
 
192
    newalign = GlobalBandToSeqAlign(gbsp);
 
193
    if (newalign != NULL)
 
194
    {    
 
195
       if (SeqLocStrand(slp1) == Seq_strand_minus) {
 
196
          sit=(SeqIntPtr)slp1->data.ptrvalue;
 
197
          sit->strand = Seq_strand_plus;
 
198
       }
 
199
       if (SeqLocStrand(slp2) == Seq_strand_minus) {
 
200
          sit=(SeqIntPtr)slp2->data.ptrvalue;
 
201
          sit->strand = Seq_strand_plus;
 
202
       } 
 
203
       AdjustOffSetsInSeqAlign (newalign, slp1, slp2);
 
204
    } 
 
205
    gbsp = DestructBandStruct(gbsp);
 
206
  } 
 
207
  else {
 
208
    ErrPostEx(SEV_WARNING, 0, 0, "Could not Create GlobalBandStruct");
 
209
  }
 
210
  return newalign;
 
211
}
 
212
 
 
213
static SeqAlignPtr BandAlignTwoSeqLocs (SeqLocPtr slp1, SeqLocPtr slp2, MashPtr msp)
 
214
{
 
215
  SeqAlignPtr         newalign = NULL;
 
216
  ValNodePtr          vnp;
 
217
 
 
218
  newalign = BandAlignTwoSeqLocsFunc (slp1, slp2, msp);
 
219
  if (newalign == NULL) {
 
220
     if ((msp->is_prot || msp->translate_prot)) {
 
221
        if (msp->reward > 1)
 
222
           msp->reward = 1;
 
223
        if (msp->penalty < -1)
 
224
           msp->penalty = -1;
 
225
     }
 
226
     else {
 
227
        if (msp->reward > 1)
 
228
           msp->reward = 1;
 
229
        if (msp->penalty < -1)
 
230
           msp->penalty = -1;
 
231
     }
 
232
     newalign = BandAlignTwoSeqLocsFunc (slp1, slp2, msp);
 
233
  } 
 
234
  if (newalign == NULL) {
 
235
     vnp = NULL;
 
236
     ValNodeAddPointer (&vnp, 0, slp1);
 
237
     ValNodeAddPointer (&vnp, 0, slp2);
 
238
     newalign = SeqLocToFastaSeqAlign (vnp); 
 
239
     ValNodeFree (vnp);
 
240
  }
 
241
  return newalign;
 
242
}
 
243
 
 
244
/********************************************************/
 
245
static Int2 number_hits (SeqAlignPtr salp)
 
246
{
 
247
  SeqAlignPtr salptmp;
 
248
  Int2 j=0;
 
249
 
 
250
  if (salp == NULL)
 
251
     return 0;
 
252
  for (salptmp=salp; salptmp!=NULL; salptmp=salptmp->next)
 
253
     j++;
 
254
  return j;
 
255
}
 
256
 
 
257
/*****************************************************
 
258
***
 
259
*** try_insert_seqalign
 
260
***    salplst = list of seqalign 
 
261
***    salp = tobe inserted
 
262
***    choice = 0 everything is kept
 
263
***             1 no overlaps
 
264
***             2 overlap in seq 1
 
265
***             3 overlap in seq 2
 
266
***             4 overlaps in both sequences of same strand
 
267
***             5 
 
268
***
 
269
*******************************************************/
 
270
static Boolean precede (Int4 pos1, Int4 pos2, Boolean is_plus)
 
271
{
 
272
  if (is_plus)
 
273
     return (Boolean) (pos1 < pos2);
 
274
  return (Boolean) (pos1 > pos2);
 
275
}
 
276
 
 
277
static SeqAlignPtr try_insert_seqalign (SeqAlignPtr salplst, SeqAlignPtr salp, Uint1 choice)
 
278
{
 
279
  SeqAlignPtr  tmp, tmp2;
 
280
  Int4         start, stop, start2, stop2,
 
281
               starti, stopi, 
 
282
               starti2, stopi2,
 
283
               startnext,
 
284
               startnext2, stopnext2;
 
285
  Boolean      goOn,
 
286
               st1, st2;
 
287
 
 
288
  if (salp == NULL)
 
289
     return NULL;
 
290
  if (salplst == NULL) {
 
291
     salplst = SeqAlignDup (salp);
 
292
     return salplst;
 
293
  }
 
294
  starti= SeqAlignStart (salp, 0);
 
295
  stopi = SeqAlignStop (salp, 0);
 
296
  starti2= SeqAlignStart (salp, 1);
 
297
  stopi2 = SeqAlignStop (salp, 1);
 
298
  st1 = (Boolean) !(SeqAlignStrand(salp,0) == Seq_strand_minus);
 
299
  st2 = (Boolean) !(SeqAlignStrand(salp,1) == Seq_strand_minus);
 
300
 
 
301
  tmp = salplst;
 
302
  start = SeqAlignStart (tmp, 0);
 
303
  stop = SeqAlignStop (tmp, 0);
 
304
  start2= SeqAlignStart (tmp, 1);
 
305
  stop2 = SeqAlignStop (tmp, 1);
 
306
  if (choice == 1)
 
307
     goOn =(Boolean) (precede(stopi, start, st1) && precede(stopi2, start2, st2));
 
308
  else if (choice == 2)
 
309
     goOn = (Boolean) (precede(stopi, start, st1) && precede(starti2, start2, st2) && precede(stopi2, stop2, st2));
 
310
  else if (choice == 3)
 
311
     goOn = (Boolean) (precede(stopi2, start2, st2) && precede(starti, start, st1) && precede(stopi, stop, st1));
 
312
  else if (choice == 4)
 
313
     goOn = (Boolean) (precede(starti, start, st1) && precede(stopi, stop, st1) && precede(starti2, start2, st2) && precede(stopi2, stop2, st2));
 
314
  else if (choice == 5)
 
315
     goOn = (Boolean) (precede(starti, start, st1));
 
316
  else
 
317
     return salplst;
 
318
  if (goOn)
 
319
  {
 
320
     tmp2 = SeqAlignDup (salp);
 
321
     tmp2->next = salplst;
 
322
     salplst = tmp2;
 
323
  }
 
324
  else {
 
325
     while (tmp != NULL)
 
326
     {
 
327
        if (choice == 1)
 
328
           goOn = (Boolean) (precede(stop, starti, st1) && precede(stop2, starti2, st2));
 
329
        else if (choice == 2)
 
330
           goOn = (Boolean) (precede(stop, starti, st1) && precede(start2, starti2, st2) && precede(stop2, stopi2, st2));
 
331
        else if (choice == 3)
 
332
           goOn = (Boolean) (precede(stop2, starti2, st2) && precede(start, starti, st1) && precede(stop, stopi, st1));
 
333
        else if (choice == 4)
 
334
           goOn = (Boolean) (precede(start, starti, st1) && precede(stop, stopi, st1) && precede(start2, starti2, st2) && precede(stop2, stopi2, st2));
 
335
        else 
 
336
           goOn=(Boolean)(precede(start, starti, st1));
 
337
        if (goOn)
 
338
        {
 
339
           if (tmp->next == NULL) {
 
340
              tmp2 = SeqAlignDup (salp);
 
341
              tmp->next = tmp2;
 
342
              break;
 
343
           }
 
344
           else {
 
345
              startnext = SeqAlignStart (tmp->next, 0);
 
346
              startnext2= SeqAlignStart (tmp->next, 1);
 
347
              stopnext2 = SeqAlignStop (tmp->next, 1);
 
348
              if (choice == 5) {
 
349
                 goOn=(Boolean)(precede(starti, startnext, st1));
 
350
              }
 
351
              else {
 
352
                 goOn=(Boolean)(precede(stopi, startnext, st1)  && precede(starti2, startnext2, st2) && precede(stopi2, stopnext2, st2)) ;
 
353
              }
 
354
              if (goOn)
 
355
              {
 
356
                 tmp2 = SeqAlignDup (salp);
 
357
                 tmp2->next = tmp->next;
 
358
                 tmp->next = tmp2;
 
359
                 break;
 
360
              }
 
361
           }
 
362
        }
 
363
        tmp = tmp->next;
 
364
        if (tmp!=NULL) {
 
365
           start = SeqAlignStart (tmp, 0);
 
366
           stop = SeqAlignStop (tmp, 0);
 
367
           start2= SeqAlignStart (tmp, 1);
 
368
           stop2 = SeqAlignStop (tmp, 1);
 
369
        }
 
370
     }
 
371
  }
 
372
  return salplst;
 
373
}
 
374
 
 
375
static SeqAlignPtr SortBlastHits (SeqAlignPtr salp, Int4 threshold, Uint1 choice)
 
376
{
 
377
  SeqAnnotPtr sap, sap2;
 
378
  SeqAlignPtr salpdup,
 
379
              salptmp,
 
380
              head = NULL,
 
381
              select = NULL,
 
382
              tmp = NULL, pre=NULL, preselect=NULL;
 
383
  Nlm_FloatHi bit_score;
 
384
  Nlm_FloatHi evalue;
 
385
  FloatLo     gap_length_ratio,
 
386
              gap_length_ratio1=1.00;
 
387
  Int4        score, 
 
388
              number,
 
389
              totalvalmax = INT4_MAX, 
 
390
              valmax;
 
391
  Int4        gap_count = 0, 
 
392
              gap_count1= 0;
 
393
  Boolean     goOn = TRUE;
 
394
  Boolean     ok;
 
395
 
 
396
  if (salp == NULL || choice < 0)
 
397
     return salp;
 
398
  sap = SeqAnnotForSeqAlign (salp);
 
399
  sap2 = (SeqAnnotPtr) AsnIoMemCopy ((Pointer) sap, (AsnReadFunc) SeqAnnotAsnRead, (AsnWriteFunc) SeqAnnotAsnWrite);
 
400
  if (sap2!=NULL) {
 
401
     salpdup = (SeqAlignPtr)sap2->data;
 
402
     sap2->data = NULL;
 
403
     SeqAnnotFree (sap2);
 
404
  }
 
405
  while (goOn) 
 
406
  { 
 
407
     valmax = 0;
 
408
     pre = select = preselect = NULL;
 
409
     for (salptmp=salpdup; salptmp!=NULL; salptmp=salptmp->next) 
 
410
     {
 
411
        GetScoreAndEvalue (salptmp, &score, &bit_score, &evalue, &number);
 
412
        gap_count = SeqAlignGapCount (salptmp);
 
413
        gap_length_ratio = (FloatLo)gap_count/SeqAlignLength(salptmp);
 
414
        if (score > valmax || (gap_count1>0 && gap_count == 0 && gap_length_ratio1 > 0.05 && score >= threshold))
 
415
        {
 
416
           valmax = score;
 
417
           select = salptmp;
 
418
           preselect = pre;
 
419
           if (head==NULL) {
 
420
              gap_count1 = gap_count; 
 
421
              gap_length_ratio1 = gap_length_ratio;
 
422
           }
 
423
        }
 
424
        pre = salptmp;
 
425
     }
 
426
     if (valmax < threshold && head == NULL)
 
427
     {
 
428
        threshold = 20;
 
429
     }
 
430
     if (valmax >= threshold && select != NULL) 
 
431
     {
 
432
        if (head!=NULL)
 
433
           ok=(Boolean)((SeqAlignStrand(head,0) == SeqAlignStrand(select,0))
 
434
              && (SeqAlignStrand(head,1) == SeqAlignStrand(select,1)));
 
435
        else ok=TRUE;
 
436
        if (ok) {
 
437
           head = try_insert_seqalign (head, select, choice);
 
438
        }
 
439
        tmp=select;
 
440
        if (preselect==NULL) {
 
441
           salpdup = select->next;
 
442
        } else {
 
443
           preselect->next = select->next; 
 
444
        }   
 
445
        tmp->next = NULL;
 
446
        SeqAlignSetFree (tmp);
 
447
        totalvalmax = valmax;
 
448
     }
 
449
     else 
 
450
        goOn = FALSE;
 
451
  }
 
452
  if (salpdup!=NULL) 
 
453
     SeqAlignSetFree (salpdup);
 
454
  return head;
 
455
}
 
456
 
 
457
static void check_strand_inSeqalign (SeqAlignPtr salp, Uint1 strand, Int2 index)
 
458
{
 
459
  DenseSegPtr dsp;
 
460
  Int2        j; 
 
461
  Uint1Ptr    strandp;
 
462
 
 
463
  if (salp!=NULL) {
 
464
     dsp = (DenseSegPtr) salp->segs;
 
465
     if (dsp->strands != NULL) {
 
466
        strandp = dsp->strands;
 
467
        strandp += index;
 
468
        for (j=0; j<dsp->numseg; j++) {
 
469
           if (*strandp != strand)
 
470
              *strandp = strand;
 
471
           strandp += dsp->dim;
 
472
        }
 
473
     }   
 
474
  }
 
475
}
 
476
 
 
477
static SeqAlignPtr SeqAlignList2PosStrand (SeqAlignPtr salp) 
 
478
{
 
479
  SeqAlignPtr salptmp;
 
480
  Int4Ptr     lenp;
 
481
  Int4Ptr     startp;
 
482
  Uint1Ptr    strandp;
 
483
  Int4        numseg;
 
484
  Int2        dim;
 
485
  Int4        j, k, n, tmp;
 
486
 
 
487
 
 
488
for (salptmp=salp; salptmp!=NULL; salptmp=salptmp->next)
 
489
{
 
490
  if ((SeqAlignStrand(salptmp,0) == Seq_strand_minus) 
 
491
   && (SeqAlignStrand(salptmp,1) != Seq_strand_minus))
 
492
  {
 
493
     if (salptmp->segtype == 2)
 
494
     {
 
495
         DenseSegPtr dsp;
 
496
        dsp = (DenseSegPtr) salptmp->segs;
 
497
        strandp = dsp->strands;
 
498
        numseg = dsp->numseg;
 
499
        dim = dsp->dim;
 
500
        lenp = dsp->lens;
 
501
        startp = dsp->starts;
 
502
        for (j=0; j < numseg*dim && strandp!=NULL; j++, strandp++)
 
503
            {
 
504
                if (*strandp == Seq_strand_minus)
 
505
                    *strandp = Seq_strand_plus;
 
506
                else if (*strandp != Seq_strand_minus)
 
507
                    *strandp = Seq_strand_minus;
 
508
            }
 
509
        for (j=0, k=numseg-1; j<numseg/2; j++, k--) {
 
510
            tmp=lenp[j];
 
511
            lenp[j]=lenp[k];
 
512
            lenp[k]=tmp;
 
513
        }
 
514
        for (j=0, k=(dim*numseg-dim); j<(dim*numseg-1)/2; j+=dim, k-=dim) {
 
515
            for (n=0; n<dim; n++) {
 
516
                tmp=startp[j+n];
 
517
                startp[j+n]=startp[k+n];
 
518
                startp[k+n]=tmp;
 
519
            }
 
520
        }
 
521
     } 
 
522
  }
 
523
}
 
524
return salp;
 
525
}
 
526
 
 
527
static SeqAlignPtr BlastTwoSeqLocs (SeqLocPtr slp1, SeqLocPtr slp2, Boolean is_prot, MashPtr msp)
 
528
{
 
529
  BLAST_OptionsBlkPtr options;
 
530
  SeqAlignPtr seqalign;
 
531
  SeqIntPtr   sit;
 
532
  Uint1       strand1, strand2;
 
533
  Boolean     delete_mash = FALSE;
 
534
 
 
535
  if (slp1 == NULL || slp2 == NULL)
 
536
     return NULL;
 
537
  if (msp == NULL) {
 
538
     msp = MashNew (is_prot);
 
539
     delete_mash = TRUE;
 
540
  }
 
541
  if (msp == NULL)
 
542
     return NULL;
 
543
  options = BLASTOptionNew((is_prot == FALSE) ? "blastn":"blastp",(Boolean) msp->use_gapped_blast);
 
544
  if (msp!=NULL) {
 
545
     options->gap_open= msp->gap_open;
 
546
     options->gap_extend= msp->gap_extend;
 
547
     if(is_prot) {
 
548
        if(msp->matrixname!=NULL) 
 
549
           options->matrix=StringSave(msp->matrixname);
 
550
     }
 
551
     else {
 
552
        options->reward = msp->reward;
 
553
        options->penalty= msp->penalty;
 
554
        options->wordsize = msp->wordsize;
 
555
     }
 
556
/*
 
557
     options->filter = 0;                 *//** msp->filter; **/
 
558
     options->gap_x_dropoff= msp->gap_x_dropoff;
 
559
     options->gap_x_dropoff_final= msp->gap_x_dropoff_final;
 
560
  }
 
561
  if (is_prot)
 
562
     options = BLASTOptionValidate (options, "blastp");
 
563
  else
 
564
     options = BLASTOptionValidate (options, "blastn");
 
565
  if (options == NULL)
 
566
     return NULL;
 
567
  if (!is_prot) {
 
568
     strand1=SeqLocStrand(slp1);
 
569
     strand2=SeqLocStrand(slp2);
 
570
     sit=(SeqIntPtr)slp1->data.ptrvalue;
 
571
     sit->strand = Seq_strand_both;
 
572
     sit=(SeqIntPtr)slp2->data.ptrvalue;
 
573
     sit->strand = Seq_strand_both;
 
574
  }
 
575
  seqalign = BlastTwoSequencesByLoc (slp1, slp2, NULL, options);
 
576
  if (options->wordsize==0) {
 
577
     while (seqalign==NULL && (is_prot==FALSE && options->wordsize>8)) {
 
578
        options->wordsize--;
 
579
        seqalign = BlastTwoSequencesByLoc (slp1, slp2, NULL, options);
 
580
     }
 
581
  }
 
582
  options=BLASTOptionDelete(options);
 
583
  if (delete_mash)
 
584
     MemFree (msp);
 
585
  if (!is_prot) {
 
586
     sit=(SeqIntPtr)slp1->data.ptrvalue;
 
587
     sit->strand = strand1;
 
588
     sit=(SeqIntPtr)slp2->data.ptrvalue;
 
589
     sit->strand = strand2;
 
590
  }
 
591
  seqalign = SeqAlignList2PosStrand (seqalign);
 
592
  return seqalign;
 
593
}
 
594
 
 
595
 
 
596
static Uint1 find_bestframe (SeqLocPtr slp, Int2 genCode)
 
597
{
 
598
  ByteStorePtr  bs;
 
599
  Int4          len;
 
600
  Int4          max = 0;
 
601
  Int4          bslen;
 
602
  Uint1         bestframe = 0, 
 
603
                frame;
 
604
  CharPtr       str;
 
605
 
 
606
  /* Only works for genCode == ncbieaa */
 
607
  for (frame = 1; frame <= 3; frame ++) 
 
608
  {
 
609
     bs = cds_to_pept (slp, frame, (Int2) Seq_code_ncbieaa, TRUE);
 
610
     if (bs!=NULL) {
 
611
        str = (CharPtr) BSMerge (bs, NULL);
 
612
        bslen=BSLen(bs);
 
613
        BSFree (bs);
 
614
        for (len = 0; len<bslen; len++)
 
615
           if (str[len]=='*')
 
616
              break;
 
617
        MemFree (str);
 
618
        if (len > max) {
 
619
           max = len;
 
620
           bestframe = frame;
 
621
        }
 
622
     }
 
623
  }
 
624
  return bestframe;
 
625
}
 
626
 
 
627
static SeqLocPtr TranslateSeqLoc (SeqLocPtr slp, Int2 genCode, Uint1 *frame)
 
628
{
 
629
  BioseqPtr     bsp = NULL;
 
630
  ByteStorePtr  bs;
 
631
  SeqLocPtr     slpp;
 
632
  Int4          length,bslen;
 
633
  
 
634
  slp = seqloc2fuzzloc (slp, TRUE, TRUE);
 
635
  *frame = find_bestframe (slp, genCode);
 
636
  bs = cds_to_pept (slp, *frame, (Int2) Seq_code_ncbieaa, TRUE);
 
637
 
 
638
  if(genCode != (Int2) Seq_code_ncbieaa ) {
 
639
      bslen=BSLen(bs);
 
640
      bs=BSConvertSeq(bs, (Uint1) genCode,(Uint1)Seq_code_ncbieaa,(Int4) bslen);
 
641
  }
 
642
  bsp = BioseqNew ();
 
643
  if (bsp != NULL) {
 
644
     bsp->repr = Seq_repr_raw;
 
645
     bsp->mol = Seq_mol_aa;
 
646
     bsp->seq_data_type = (Uint1)genCode;
 
647
     bsp->seq_data = bs;
 
648
     bsp->length = BSLen (bs);
 
649
     bsp->id = MakeNewProteinSeqId (slp, NULL);
 
650
     if (*frame==2)
 
651
        length = (Int4)((SeqLocLen(slp)-1) / (Int4)3);
 
652
     else if (*frame==3) 
 
653
        length = (Int4)((SeqLocLen(slp)-2) / (Int4)3);
 
654
     else
 
655
        length = (Int4)(SeqLocLen(slp) / (Int4)3);
 
656
     slpp = SeqLocIntNew (0, length-1, Seq_strand_plus, bsp->id);
 
657
  }
 
658
  return slpp;
 
659
}
 
660
 
 
661
 
 
662
 
 
663
 
 
664
/*************************************************************
 
665
*** GlobalAlignTwoSeqLocs
 
666
***   if short sequences (length < 150)
 
667
***         aligns using BandAlign only
 
668
***   else 
 
669
***         Blast the 2 seqlocs
 
670
***         if find no hit:
 
671
***            aligns using BandAlign
 
672
***         else if finds 1 hit (NOW: the longest)
 
673
***            if alignment reaches ends (3' or 5')
 
674
***               extend seqalign with the missing ends
 
675
***            else 
 
676
***               aligns unaligned seqlocs using BandAlign
 
677
***               merge seqaligns
 
678
***
 
679
*** !!!SelectBestHit : select 1 hit
 
680
***                    should select several if any
 
681
***
 
682
**************************************************************/
 
683
static SeqAlignPtr AlignExtreme5 (SeqAlignPtr salp, MashPtr msp, Int4 slpstart1, Int4 start1, Int4 slpstart2, Int4 slpstop2, Int4 start2, SeqIdPtr sip1, SeqIdPtr sip2, Uint1 strand1, Uint1 strand2)
 
684
{
 
685
  SeqAlignPtr salp2;
 
686
  SeqLocPtr   newslp1, newslp2;
 
687
  Boolean     touch_end5;
 
688
  
 
689
  if (strand2 == Seq_strand_minus) 
 
690
  {
 
691
         touch_end5 = (Boolean)((slpstart1==start1) || (slpstop2==start2));
 
692
         if (touch_end5) 
 
693
         {
 
694
            salp = SeqAlignEndExtend (salp, slpstart1, slpstop2, -1, -1, start1, start2, -1, -1, strand1, strand2);
 
695
         } 
 
696
         else 
 
697
         {
 
698
            start1 -= 1; 
 
699
            start2 += 1; 
 
700
            newslp1=SeqLocIntNew (0, start1, strand1, sip1);
 
701
            newslp2=SeqLocIntNew (start2, slpstop2, strand2, sip2);
 
702
            msp->rg1_ext=(Int4)msp->gap_extend;
 
703
            msp->lg1_ext=0;
 
704
            msp->rg2_ext=(Int4)msp->gap_extend;
 
705
            msp->lg2_ext=0;
 
706
            msp->rg1_open=(Int4)msp->gap_open;
 
707
            msp->lg1_open=0;
 
708
            msp->rg2_open=(Int4)msp->gap_open;
 
709
            msp->lg2_open=0;
 
710
            salp2 = BandAlignTwoSeqLocs (newslp1, newslp2, msp);
 
711
            salp = SeqAlignMerge (salp, salp2, FALSE);
 
712
         }
 
713
  } 
 
714
  else 
 
715
  {
 
716
         touch_end5 = (Boolean)((ABS(slpstart1-start1)<3) || (ABS(slpstart2-start2)<2));
 
717
         if (touch_end5) 
 
718
         {
 
719
            salp = SeqAlignEndExtend (salp, slpstart1, slpstart2, -1, -1, start1, start2, -1, -1, strand1, strand2);
 
720
         } 
 
721
         else 
 
722
         {
 
723
            start1 -= 1; 
 
724
            start2 -= 1; 
 
725
            newslp1=SeqLocIntNew (0, start1, strand1, sip1);
 
726
            newslp2=SeqLocIntNew (0, start2, strand2, sip2);
 
727
            msp->rg1_ext=(Int4)msp->gap_extend;
 
728
            msp->lg1_ext=0;
 
729
            msp->rg2_ext=(Int4)msp->gap_extend;
 
730
            msp->lg2_ext=0;
 
731
            msp->rg1_open=(Int4)msp->gap_open;
 
732
            msp->lg1_open=0;
 
733
            msp->rg2_open=(Int4)msp->gap_open;
 
734
            msp->lg2_open=0;
 
735
            salp2 = BandAlignTwoSeqLocs (newslp1, newslp2, msp);
 
736
            salp = SeqAlignMerge (salp, salp2, FALSE);
 
737
         }
 
738
  }
 
739
  return salp;
 
740
}
 
741
 
 
742
static SeqAlignPtr AlignExtreme3 (SeqAlignPtr salp, MashPtr msp, Int4 slpstop1, Int4 stop1, Int4 slpstart2, Int4 slpstop2, Int4 stop2, SeqIdPtr sip1, SeqIdPtr sip2, Uint1 strand1, Uint1 strand2)
 
743
{
 
744
  SeqAlignPtr salp2;
 
745
  SeqLocPtr   newslp1, newslp2;
 
746
  Boolean     touch_end3;
 
747
  
 
748
  if (strand2 == Seq_strand_minus) 
 
749
  {
 
750
         touch_end3 = (Boolean)((slpstop1==stop1) || (slpstart2==stop2));
 
751
         if (touch_end3) 
 
752
         {
 
753
            salp = SeqAlignEndExtend (salp, -1, -1, slpstop1, slpstart2, -1, -1, stop1, stop2, strand1, strand2);
 
754
         } 
 
755
         else 
 
756
         {
 
757
            stop1 += 1;
 
758
            stop2 -= 1;
 
759
            newslp1=SeqLocIntNew(stop1,slpstop1, strand1, sip1);
 
760
            newslp2=SeqLocIntNew(0,stop2, strand2, sip2);
 
761
            msp->lg1_ext= (Int4)msp->gap_extend;
 
762
            msp->rg1_ext=0;
 
763
            msp->lg2_ext= (Int4) msp->gap_extend;
 
764
            msp->rg2_ext=0;
 
765
            msp->lg1_open= (Int4)msp->gap_open;
 
766
            msp->rg1_open=0;
 
767
            msp->lg2_open= (Int4)msp->gap_open;
 
768
            msp->rg2_open=0;
 
769
            salp2= BandAlignTwoSeqLocs (newslp1,newslp2, msp);
 
770
            salp = SeqAlignMerge (salp, salp2, TRUE);
 
771
         } 
 
772
  } 
 
773
  else 
 
774
  {
 
775
         touch_end3 = (Boolean)((ABS(slpstop1-stop1)<3) || (ABS(slpstop2-stop2)<3));
 
776
         if (touch_end3) 
 
777
         {
 
778
            salp = SeqAlignEndExtend (salp, -1, -1, slpstop1, slpstop2, -1, -1, stop1, stop2, strand1, strand2);
 
779
         } 
 
780
         else 
 
781
         {
 
782
            stop1 += 1;
 
783
            stop2 += 1;
 
784
            newslp1=SeqLocIntNew(stop1,slpstop1, strand1, sip1);
 
785
            newslp2=SeqLocIntNew(stop2,slpstop2, strand2, sip2);
 
786
            msp->lg1_ext=(Int4)msp->gap_extend;
 
787
            msp->rg1_ext=0;
 
788
            msp->lg2_ext=(Int4) msp->gap_extend;
 
789
            msp->rg2_ext=0;
 
790
            msp->lg1_open=(Int4)msp->gap_open;
 
791
            msp->rg1_open=0;
 
792
            msp->lg2_open=(Int4)msp->gap_open;
 
793
            msp->rg2_open=0;
 
794
            salp2 = BandAlignTwoSeqLocs (newslp1, newslp2, msp);
 
795
            salp = SeqAlignMerge (salp, salp2, TRUE);
 
796
         } 
 
797
  }
 
798
  return salp;
 
799
}
 
800
 
 
801
static SeqAlignPtr align_extrem (SeqAlignPtr salp, SeqLocPtr slp1, SeqLocPtr slp2, MashPtr msp)
 
802
{
 
803
  SeqAlignPtr tmp;
 
804
  SeqIdPtr    sip1, sip2;
 
805
  Int4        start1, start2,
 
806
              stop1, stop2;
 
807
  Int4        slpstart1, slpstart2,
 
808
              slpstop1, slpstop2;
 
809
  Uint1       strand1, strand2;
 
810
  
 
811
  sip1 = SeqLocId(slp1);
 
812
  sip2 = SeqLocId(slp2);
 
813
  strand1 = SeqAlignStrand (salp, 0);
 
814
  strand2 = SeqAlignStrand (salp, 1);
 
815
  slpstart1= SeqLocStart(slp1);
 
816
  slpstop1 = SeqLocStop(slp1);
 
817
  slpstart2= SeqLocStart(slp2);
 
818
  slpstop2 = SeqLocStop(slp2);  
 
819
  start1 = SeqAlignStart (salp, 0);
 
820
  start2 = SeqAlignStart (salp, 1);
 
821
  salp = AlignExtreme5 (salp, msp, slpstart1, start1, slpstart2, slpstop2, start2, sip1, sip2, strand1, strand2);
 
822
  check_strand_inSeqalign (salp, strand1, 0);
 
823
  check_strand_inSeqalign (salp, strand2, 1);
 
824
 
 
825
  if (salp!=NULL)
 
826
  {
 
827
     tmp = salp;
 
828
     while (tmp->next!=NULL)
 
829
        tmp = tmp->next;
 
830
     stop1 = SeqAlignStop (tmp, 0);
 
831
     stop2 = SeqAlignStop (tmp, 1);  
 
832
     tmp = AlignExtreme3 (tmp, msp, slpstop1,  stop1,  slpstart2, slpstop2, stop2, sip1, sip2, strand1, strand2);
 
833
     check_strand_inSeqalign (tmp, strand1, 0);
 
834
     check_strand_inSeqalign (tmp, strand2, 1);
 
835
  }
 
836
  return salp;
 
837
}
 
838
 
 
839
/******************************************************
 
840
static Boolean is_intron (CharPtr str, Int4 len)
 
841
{
 
842
  if (str[0]=='G' && str[1]=='T' && str[len-2]=='A' && str[len-1]=='G')
 
843
     return TRUE;
 
844
  return FALSE;
 
845
}
 
846
*********************************************************/
 
847
static FloatHi is_donor (CharPtr str, Int4 len)
 
848
{
 
849
  FloatHi one[4]={0.35, 0.35, 0.19, 0.11};
 
850
  FloatHi two[4]={0.59, 0.13, 0.14, 0.14};
 
851
  FloatHi three[4]={0.08, 0.02, 0.82, 0.08};
 
852
  FloatHi four[4]={0.01, 0.01, 1.00, 0.01};
 
853
  FloatHi five[4]={0.01, 0.01, 0.01, 1.00};
 
854
  FloatHi six[4]={0.51, 0.03, 0.43, 0.03};
 
855
  FloatHi seven[4]={0.71, 0.08, 0.12, 0.09};
 
856
  FloatHi eight[4]={0.06, 0.05, 0.84, 0.05};
 
857
  FloatHi nine[4]={0.15, 0.16, 0.17, 0.52};
 
858
  FloatHi score =1.000;
 
859
  Int4   i;
 
860
  Int4  PNTR num=NULL;
 
861
 
 
862
  if ((num = (Int4Ptr)MemNew(len*sizeof(Int4)))==NULL) {
 
863
    return(-1);
 
864
  }
 
865
 
 
866
  for (i = 0; i <= len; i++){
 
867
    if (str[i]=='A')
 
868
      num[i] = 0;
 
869
    if (str[i]=='C')
 
870
      num[i] = 1;
 
871
    if (str[i]=='G')
 
872
      num[i] = 2;
 
873
    if (str[i]=='T')
 
874
      num[i] = 3;
 
875
  }
 
876
  score *= one[num[0]];
 
877
  score *= two[num[1]];
 
878
  score *= three[num[2]];
 
879
  score *= four[num[3]];
 
880
  score *= five[num[4]];
 
881
  score *= six[num[5]];
 
882
  score *= seven[num[6]];
 
883
  score *= eight[num[7]];
 
884
  score *= nine[num[8]];
 
885
 
 
886
  MemFree(num);
 
887
  num=NULL;
 
888
 
 
889
  return score;
 
890
}
 
891
 
 
892
static Int4 getSplicePos (CharPtr str)
 
893
{
 
894
  Int4     offset = -1;
 
895
  Int4     xcursor = 0;
 
896
  Int4     c;
 
897
  FloatHi  topscore = -FLT_MAX,
 
898
           score;
 
899
  Char     tmpstr[9];
 
900
  Int4     length;
 
901
 
 
902
  if (str == NULL)
 
903
    return -1;
 
904
  length = MIN(StringLen(str)/2-10, 10);
 
905
  while (xcursor <= length)
 
906
  {
 
907
      for (c = 0; c < 9; c++)
 
908
      {
 
909
        tmpstr[c] = str[xcursor+c];
 
910
      }
 
911
      if ((score=is_donor(tmpstr, 8)) > topscore)
 
912
      {
 
913
        topscore = score;
 
914
        offset = xcursor;
 
915
      }
 
916
      xcursor += 1;
 
917
  }
 
918
  if (topscore > 0.000010 && offset >=0)
 
919
    return offset+3;
 
920
  return -1;
 
921
}
 
922
 
 
923
 
 
924
 
 
925
static SeqAlignPtr align_btwhits (SeqAlignPtr salp, SeqIdPtr sip1, SeqIdPtr sip2, MashPtr msp)
 
926
{
 
927
  SeqAlignPtr tmp,
 
928
              newsalp,
 
929
              newsalphead, newtmp;
 
930
  SeqLocPtr   slp1, slp2;
 
931
  Int4        x1, y1, x2, y2;
 
932
  Int4        len;
 
933
  Int4        offset;
 
934
  Uint1       st1, st2;
 
935
  CharPtr     str;
 
936
  Boolean     search_intron = FALSE;
 
937
 
 
938
  if (salp == NULL) return NULL;
 
939
  if (msp != NULL) {
 
940
     search_intron = msp->splicing;
 
941
  }
 
942
  newsalphead = SeqAlignDup (salp);
 
943
  newtmp = newsalphead;
 
944
  tmp = salp ->next;
 
945
  while (tmp != NULL)
 
946
  {
 
947
     x1 = SeqAlignStop (newtmp, 0);
 
948
     y1 = SeqAlignStart (tmp, 0);
 
949
     x2 = SeqAlignStop (newtmp, 1);
 
950
     y2 = SeqAlignStart (tmp, 1);
 
951
     st1= SeqAlignStrand (newtmp, 0);
 
952
     st2= SeqAlignStrand (newtmp, 1);
 
953
 
 
954
     if (x1 + 1 == y1) {
 
955
        if (y2 == x2 + 1)
 
956
        {
 
957
           newtmp = SeqAlignMerge (newtmp, tmp, TRUE);
 
958
           tmp = tmp->next;
 
959
        }
 
960
        else if (x2 >= y2) {
 
961
           SeqAlignTrunc2 (newtmp, 0, -(abs(x2-y2+1)));
 
962
        }
 
963
        else {
 
964
           if(st1!=Seq_strand_minus) y1 -= 1; else y1+=1;
 
965
           if(st2!=Seq_strand_minus) y2 -= 1; else y2+=1;
 
966
           newtmp = SeqAlignEndExtend (newtmp, -1, -1, y1, y2, -1, -1, x1, x2, st1, st2);
 
967
           newtmp = SeqAlignMerge (newtmp, tmp, TRUE);
 
968
           tmp = tmp->next;
 
969
        }
 
970
     }
 
971
     else if (x1 >= y1) {
 
972
        SeqAlignTrunc2 (newtmp, +(abs(x1-y1+1)), 0);
 
973
     }
 
974
     else
 
975
     {
 
976
        if ((st2!=Seq_strand_minus && y2 == x2 + 1) || (st2==Seq_strand_minus && x2 == y2+1))
 
977
        {
 
978
           if (search_intron) 
 
979
           {
 
980
              str =load_seq_data(sip1, x1-5, y1+1, FALSE, &len);
 
981
              offset = getSplicePos (str);
 
982
              if (offset>= -1 && offset<len-1) 
 
983
              {
 
984
                 offset -= 6;
 
985
                 if ((offset>=-6 && offset<0) || (offset>0)) 
 
986
                 {
 
987
                    SeqAlignTrunc2 (newtmp, 0, abs(offset));
 
988
                    SeqAlignTrunc2 (tmp, abs(offset), 0);
 
989
                 }
 
990
                 x1 = SeqAlignStop (newtmp, 0);
 
991
                 y1 = SeqAlignStart (tmp, 0);
 
992
                 x2 = SeqAlignStop (newtmp, 1);
 
993
                 y2 = SeqAlignStart (tmp, 1);
 
994
                 if(st1!=Seq_strand_minus) y1 -= 1; else y1+=1;
 
995
                 if(st2!=Seq_strand_minus) y2 -= 1; else y2+=1;
 
996
                 newtmp=SeqAlignEndExtend (newtmp,-1,-1,y1, y2, -1, -1, x1, x2, st1, st2);
 
997
                 newtmp = SeqAlignMerge (newtmp, tmp, TRUE);
 
998
                 tmp = tmp->next;
 
999
              }
 
1000
              MemFree (str);
 
1001
           } 
 
1002
           else  {
 
1003
              if(st1!=Seq_strand_minus) y1 -= 1; else y1+=1; 
 
1004
              if(st2!=Seq_strand_minus) y2 -= 1; else y2+=1; 
 
1005
              newtmp=SeqAlignEndExtend (newtmp, -1, -1, y1, y2, -1, -1, x1, x2, st1, st2);
 
1006
              newtmp = SeqAlignMerge (newtmp, tmp, TRUE);
 
1007
              tmp = tmp->next;
 
1008
           }
 
1009
        }
 
1010
        else if (st2!=Seq_strand_minus && x2 >= y2) {
 
1011
            SeqAlignTrunc2 (newtmp, 0, -(abs(x2-y2+1)));
 
1012
        }
 
1013
        else if (st2==Seq_strand_minus && y2 >= x2) {
 
1014
            SeqAlignTrunc2 (tmp, abs(y2-x2+1), 0);
 
1015
        }
 
1016
        else {
 
1017
           slp1 = SeqLocIntNew (x1+1, y1-1, st1, sip1);
 
1018
           if (st2!=Seq_strand_minus)
 
1019
              slp2 = SeqLocIntNew (x2+1, y2-1, st2, sip2);
 
1020
           else 
 
1021
              slp2 = SeqLocIntNew (y2+1, x2-1, st2, sip2);
 
1022
           newsalp = BandAlignTwoSeqLocs (slp1, slp2, msp);
 
1023
           if (newsalp != NULL) 
 
1024
           {
 
1025
              newtmp = SeqAlignMerge (newtmp, newsalp, TRUE);
 
1026
              newtmp = SeqAlignMerge (newtmp, tmp, TRUE);
 
1027
           } 
 
1028
           tmp = tmp->next;     
 
1029
        } 
 
1030
     }
 
1031
  }
 
1032
  return newsalphead;
 
1033
}
 
1034
 
 
1035
/*************************************/
 
1036
static Boolean possible_child (SeqAlignPtr salp1, SeqAlignPtr salp2, Uint1 choice)
 
1037
{
 
1038
  Int4         start, stop, start2, stop2,
 
1039
               starti, stopi, 
 
1040
               starti2, stopi2;
 
1041
  Boolean      goOn,
 
1042
               st1, st2;
 
1043
 
 
1044
  starti= SeqAlignStart (salp1, 0);
 
1045
  stopi = SeqAlignStop (salp1, 0);
 
1046
  starti2= SeqAlignStart (salp1, 1);
 
1047
  stopi2 = SeqAlignStop (salp1, 1);
 
1048
  st1 = (Boolean) !(SeqAlignStrand(salp1,0) == Seq_strand_minus);
 
1049
  st2 = (Boolean) !(SeqAlignStrand(salp1,1) == Seq_strand_minus);
 
1050
 
 
1051
  start = SeqAlignStart (salp2, 0);
 
1052
  stop = SeqAlignStop (salp2, 0);
 
1053
  start2= SeqAlignStart (salp2, 1);
 
1054
  stop2 = SeqAlignStop (salp2, 1);
 
1055
  if (choice == 1)
 
1056
     goOn = (Boolean) (precede(stopi, start, st1) && precede(stopi2, start2, st2));
 
1057
  else if (choice == 2)
 
1058
     goOn = (Boolean) (precede(stopi, start, st1) && precede(starti2, start2, st2) && precede(stopi2, stop2, st2));
 
1059
  else if (choice == 3)
 
1060
     goOn = (Boolean) (precede(stopi2, start2, st2) && precede(starti, start, st1) && precede(stopi, stop, st1));
 
1061
  else if (choice == 4)
 
1062
     goOn = (Boolean) (precede(starti, start, st1) && precede(stopi, stop, st1) && precede(starti2, start2, st2) && precede(stopi2, stop2, st2));
 
1063
  else if (choice == 5)
 
1064
     goOn = (Boolean) (precede(starti, start, st1));
 
1065
  else
 
1066
     goOn = TRUE;
 
1067
  return goOn;
 
1068
}
 
1069
 
 
1070
static NodeHitPtr CreateGraph (SeqAlignPtr salp, Uint1 choice)
 
1071
{
 
1072
  NodeHitPtr    headnode, 
 
1073
                node1, node2,
 
1074
                newnode,
 
1075
                child;
 
1076
  SeqAlignPtr   tmp1, tmp2;
 
1077
  Nlm_FloatHi   bit_score;
 
1078
  Nlm_FloatHi   evalue;
 
1079
  Int4          score, number;
 
1080
  Int2          j;
 
1081
  
 
1082
  headnode = (NodeHitPtr) MemNew (sizeof(NodeHit));
 
1083
  headnode->index = 1;
 
1084
  headnode->salp = salp;
 
1085
  headnode->child = NULL;
 
1086
  headnode->next = NULL;
 
1087
  node1 = headnode;
 
1088
  for (tmp1=salp->next, j=2; tmp1!=NULL; tmp1=tmp1->next, j++) 
 
1089
  {
 
1090
     newnode = (NodeHitPtr) MemNew (sizeof(NodeHit));
 
1091
     newnode->index = j;
 
1092
     newnode->open = TRUE;
 
1093
     newnode->salp = tmp1;
 
1094
     newnode->child = NULL;
 
1095
     newnode->next = NULL;
 
1096
     GetScoreAndEvalue (tmp1, &score, &bit_score, &evalue, &number);
 
1097
     newnode->evalue = evalue;
 
1098
     newnode->score = score;
 
1099
     newnode->bit_score = bit_score;
 
1100
     node1->next = newnode;
 
1101
     node1 = node1->next;
 
1102
  }
 
1103
  node1 = headnode;
 
1104
  for (tmp1=salp; tmp1!=NULL; tmp1=tmp1->next) 
 
1105
  {
 
1106
     child = NULL;
 
1107
     node2 = headnode;
 
1108
     for (tmp2=salp; tmp2!=NULL; tmp2=tmp2->next) 
 
1109
     {
 
1110
        if (possible_child (tmp1, tmp2, choice)) {
 
1111
           newnode = (NodeHitPtr) MemNew (sizeof(NodeHit));
 
1112
           newnode->index = -1;
 
1113
           newnode->open = TRUE;
 
1114
           newnode->salp = NULL;
 
1115
           newnode->child = node2;
 
1116
           newnode->next = NULL;
 
1117
           if (node1->child == NULL) {
 
1118
              node1->child = newnode;
 
1119
              child = node1->child;
 
1120
           } else {
 
1121
              child->next = newnode;
 
1122
              child = child->next;
 
1123
           }
 
1124
        }
 
1125
        node2 = node2->next;
 
1126
     }
 
1127
     node1 = node1->next;
 
1128
  }
 
1129
  return headnode;
 
1130
}
 
1131
 
 
1132
static NodeHitPtr SimplifyGraph (NodeHitPtr headnode)
 
1133
{
 
1134
  NodeHitPtr node1, node2, node3, node4, node5;
 
1135
  Int2       gdchild;
 
1136
 
 
1137
  for (node1 = headnode; node1!=NULL; node1 = node1->next)
 
1138
  {
 
1139
     for (node2=node1->child; node2!=NULL; node2 = node2->next)
 
1140
     {
 
1141
        node3=node2->child;
 
1142
        if (node3) {
 
1143
           for (node4=node3->child; node4!=NULL; node4=node4->next)
 
1144
           {
 
1145
              gdchild = node4->child->index;
 
1146
              for (node5=node1->child; node5!=NULL; node5 = node5->next)
 
1147
              {
 
1148
                 if (node5->child->index == gdchild)
 
1149
                 {
 
1150
                    node5->open = FALSE;
 
1151
                    break;
 
1152
                 }
 
1153
              }
 
1154
           }
 
1155
        }
 
1156
     }
 
1157
  }
 
1158
  return headnode;
 
1159
}
 
1160
 
 
1161
 
 
1162
static SeqAlignPtr link_solution (ValNodePtr vnp, NodeHitPtr head, Int2 total)
 
1163
{
 
1164
  NodeHitPtr  nodetmp;
 
1165
  SeqAlignPtr headsalp = NULL,
 
1166
              tmpsalp,
 
1167
              newsalp, salp;
 
1168
  Int2Ptr     tab;
 
1169
  Int2        j;
 
1170
  Boolean     first = TRUE;
 
1171
 
 
1172
  if (vnp==NULL)
 
1173
     return NULL;
 
1174
  tab = (Int2Ptr) vnp->data.ptrvalue;
 
1175
  if (tab==NULL)
 
1176
     return NULL;
 
1177
  for (j=0; j<total; j++) 
 
1178
  {
 
1179
     if (tab[j] > -1) 
 
1180
     {
 
1181
        for(nodetmp=head; nodetmp!=NULL; nodetmp=nodetmp->next) {
 
1182
           if (nodetmp->index == tab[j])
 
1183
              break;
 
1184
        }
 
1185
        if (nodetmp!=NULL && nodetmp->index == tab[j]) {
 
1186
           salp=nodetmp->salp;
 
1187
           if (salp!=NULL) {
 
1188
              newsalp = SeqAlignDup (salp);
 
1189
              if (headsalp==NULL) {
 
1190
                 headsalp = newsalp;
 
1191
                 tmpsalp = headsalp;
 
1192
              } else {
 
1193
                 tmpsalp->next = newsalp;
 
1194
              }
 
1195
              tmpsalp = newsalp;
 
1196
           }
 
1197
        }
 
1198
     }
 
1199
     else break;
 
1200
  }
 
1201
  return headsalp;
 
1202
}
 
1203
 
 
1204
static ValNodePtr find_maxsolution (ValNodePtr vnp, NodeHitPtr head, Int2 total, float *maxscore, float *minintron, Int4 *alignlens)
 
1205
{
 
1206
  ValNodePtr  tmp;
 
1207
  NodeHitPtr  nodetmp;
 
1208
  SeqAlignPtr salp;
 
1209
  Int2Ptr     tab;
 
1210
  float       sum;
 
1211
  Int4        start, start1, stop;
 
1212
  Int4        alignlen;
 
1213
  Int4        bestlen = 0;
 
1214
  float       bestscore = -100.00;
 
1215
  float       intronlg;
 
1216
  float       bestintron = FLT_MAX;
 
1217
  Int2        j;
 
1218
  Boolean     first=TRUE;
 
1219
 
 
1220
  for (tmp=vnp; tmp!=NULL; tmp=tmp->next)
 
1221
  {
 
1222
     tab = (Int2Ptr) tmp->data.ptrvalue;
 
1223
     sum = 0;
 
1224
     intronlg = (float)0.00;
 
1225
     alignlen = 0;
 
1226
     for (j=0; j<total; j++) 
 
1227
     {
 
1228
      if (tab[j] > -1) 
 
1229
      {
 
1230
         for(nodetmp=head; nodetmp!=NULL; nodetmp=nodetmp->next) {
 
1231
            if (nodetmp->index == tab[j])
 
1232
               break;
 
1233
         }
 
1234
          if (nodetmp!=NULL && nodetmp->index == tab[j]) {
 
1235
            salp=nodetmp->salp;
 
1236
            if (salp!=NULL) {
 
1237
               sum += nodetmp->score;
 
1238
               start = SeqAlignStart(salp,0);
 
1239
               if (first) {
 
1240
                  first = FALSE;
 
1241
                  start1 = start;
 
1242
               }
 
1243
               stop = SeqAlignStop(salp,0);
 
1244
               alignlen += abs(stop - start);
 
1245
/***
 
1246
WriteLog ("%ld..%ld %ld %ld ", start, stop, alignlen, (long)sum);
 
1247
***/
 
1248
            }
 
1249
         }
 
1250
      }
 
1251
      else {
 
1252
         intronlg = (float)(abs(stop - start1) - alignlen)/(float)(j-1);
 
1253
         if (sum > bestscore) {
 
1254
            bestscore = sum;
 
1255
         }
 
1256
         if (intronlg < bestintron) {
 
1257
            bestintron = intronlg;
 
1258
         }
 
1259
         if (alignlen > bestlen) {
 
1260
            bestlen = alignlen;
 
1261
         }
 
1262
/***
 
1263
WriteLog ("= %ld > %ld, %f > %f, %ld > %ld\n",(long)sum, (long)bestscore, intronlg, bestintron, alignlen, bestlen);
 
1264
***/
 
1265
         break;
 
1266
      }
 
1267
     }
 
1268
  }
 
1269
  *maxscore = bestscore;
 
1270
  *minintron = bestintron;
 
1271
  *alignlens = bestlen;
 
1272
  return vnp;
 
1273
}
 
1274
 
 
1275
static ValNodePtr get_solutions (ValNodePtr vnp, NodeHitPtr head, Int2 total, Int4 totalens)
 
1276
{
 
1277
  ValNodePtr  tmp, 
 
1278
              bestvnp;
 
1279
  NodeHitPtr  nodetmp;
 
1280
  SeqAlignPtr salp;
 
1281
  Int2Ptr     tab;
 
1282
  float       sum;
 
1283
  Int4        start, start1, stop;
 
1284
  float       maxscore;
 
1285
  Int4        maxalignlens, alignlens;
 
1286
  Int2        index, 
 
1287
              j;
 
1288
  float       intronlg;
 
1289
  float       minintron;
 
1290
  float       bestratio=0.00;
 
1291
  Int4        bestratio1=0;
 
1292
  float       x, y, z;
 
1293
  Boolean     first=TRUE;
 
1294
 
 
1295
  find_maxsolution (vnp, head, total, &maxscore, &minintron, &maxalignlens);
 
1296
  bestvnp = NULL;
 
1297
  index = 1;
 
1298
  for (tmp=vnp; tmp!=NULL; tmp=tmp->next)
 
1299
  {
 
1300
    if(tmp->choice==0) 
 
1301
    {
 
1302
     tab = (Int2Ptr) tmp->data.ptrvalue;
 
1303
     sum = 0;
 
1304
     intronlg = (float)0.00;
 
1305
     alignlens = 0;
 
1306
     first = TRUE;
 
1307
     for (j=0; j<total; j++) 
 
1308
     {
 
1309
      if (tab[j] > -1) 
 
1310
      {
 
1311
         for(nodetmp=head; nodetmp!=NULL; nodetmp=nodetmp->next) {
 
1312
            if (nodetmp->index == tab[j])
 
1313
               break;
 
1314
         }
 
1315
         if (nodetmp!=NULL && nodetmp->index == tab[j]) {
 
1316
            salp=nodetmp->salp;
 
1317
            if (salp!=NULL) {
 
1318
               sum += nodetmp->score;
 
1319
               start = SeqAlignStart(salp,0);
 
1320
               if (first) {
 
1321
                  first = FALSE;
 
1322
                  start1 = start;
 
1323
               }
 
1324
               stop = SeqAlignStop(salp,0);
 
1325
               alignlens += abs(stop - start);
 
1326
            }
 
1327
         }
 
1328
      }
 
1329
      else {
 
1330
         intronlg = (float)(abs(stop - start1) - alignlens)/(float)(j-1);
 
1331
         x = (float)sum / (float)maxscore;
 
1332
         y = (float)intronlg / (float)minintron; 
 
1333
         z = (float)alignlens / (float)maxalignlens;
 
1334
         if ((Int4)(1000.0*x*z) > bestratio1 && (float)(x*z/y) > (float)bestratio) {
 
1335
/****
 
1336
WriteLog ("FFF %ld %ld =%f / %f %d = %f/  %d %d  %f    %f   %f %ld  ++ %d %d %d %d  %f  %f\n", (long)sum, (long)maxscore, x, intronlg, j-1, y, alignlens, maxalignlens, z, (float)((x*z)/y), (float)bestratio, (long)bestratio1, abs(stop - start1), alignlens, start1, stop, (float)(abs(stop - start1) - alignlens), (float)minintron );
 
1337
*****/
 
1338
            bestratio1 = (Int4) 1000*(Int4)(x*z);
 
1339
            bestratio = (float)(x*z/y);
 
1340
            bestvnp = tmp;
 
1341
         }
 
1342
         index++;
 
1343
         break;
 
1344
      }
 
1345
     }
 
1346
    }
 
1347
  }
 
1348
  if (bestvnp!=NULL)
 
1349
     bestvnp->choice = 5;
 
1350
  return bestvnp;
 
1351
}
 
1352
 
 
1353
static ValNodePtr traverseGraph (NodeHitPtr head, NodeHitPtr node, Int2 total, Int2 level, Int2Ptr tab, ValNodePtr vnp, Boolean * change)       
 
1354
{
 
1355
  NodeHitPtr  child,
 
1356
              nhtmp;
 
1357
  Int2Ptr     tab2;
 
1358
  Int2        j;
 
1359
 
 
1360
  child = node->child;
 
1361
  tab[level] = node->index;
 
1362
  *change = TRUE;
 
1363
 
 
1364
  while (child!=NULL) 
 
1365
  {
 
1366
     nhtmp = child->child;
 
1367
     if (child->open && nhtmp!=NULL) {
 
1368
        vnp = traverseGraph(head, nhtmp, total, level+1, tab, vnp, change);
 
1369
     }
 
1370
     child = child->next;
 
1371
  }
 
1372
  if (level > 0 && *change) 
 
1373
  {
 
1374
     tab2 = (Int2Ptr)MemNew ((size_t)((total+2)*sizeof(Int2)));
 
1375
     for (j=0; j<=level; j++) {
 
1376
        tab2[j] = tab[j];
 
1377
     }
 
1378
     for (; j<total; j++)
 
1379
        tab2 [j] = -1;
 
1380
     ValNodeAddPointer (&vnp, 0, (Pointer)tab2);
 
1381
     tab[level] = -1;
 
1382
     *change = FALSE;
 
1383
  }
 
1384
  return vnp;
 
1385
}
 
1386
 
 
1387
static SeqAlignPtr seqalign_reverse_sorting (SeqAlignPtr salp)
 
1388
{
 
1389
  SeqAlignPtr salptmp, tmp2,next,
 
1390
              tail,
 
1391
              salpnew=NULL;
 
1392
 
 
1393
  salptmp = salp;
 
1394
  while (salptmp->next)
 
1395
  {
 
1396
     tmp2=salptmp; 
 
1397
     next = tmp2->next;
 
1398
     while (next->next) {
 
1399
        tmp2=tmp2->next;
 
1400
        next=next->next;
 
1401
     }
 
1402
     if (tmp2->next!=NULL) {
 
1403
        if (salpnew==NULL) {
 
1404
           salpnew = tmp2->next;
 
1405
           tail = salpnew;
 
1406
        } else  {
 
1407
           tail->next = tmp2->next;
 
1408
           tail = tail->next;
 
1409
        }
 
1410
        tmp2->next = NULL;
 
1411
     }
 
1412
  }
 
1413
  if (salptmp)
 
1414
  {
 
1415
     tail->next = salptmp;
 
1416
  }
 
1417
  return salpnew;
 
1418
}
 
1419
 
 
1420
 
 
1421
/*************************************/
 
1422
static SeqAlignPtr BlastBandAlignTwoSeqLocs (SeqLocPtr slp1, SeqLocPtr slp2, SeqAlignPtr salp, Boolean is_prot, MashPtr msp)
 
1423
{
 
1424
  SeqAlignPtr newsalp = NULL;
 
1425
  Uint1       strand;
 
1426
  Boolean     delete_mash = FALSE;
 
1427
 
 
1428
  ValNodePtr  vnp;
 
1429
  NodeHitPtr  headnode, node;
 
1430
  ValNodePtr  bestsol;
 
1431
  Int2Ptr     tab;
 
1432
  Int2        j, total;
 
1433
  Boolean     change;
 
1434
 
 
1435
  if (salp != NULL) {
 
1436
     if (msp == NULL) {
 
1437
        msp = MashNew (FALSE);
 
1438
        delete_mash = TRUE;
 
1439
     }
 
1440
     if (msp == NULL)
 
1441
        return NULL;
 
1442
     if (number_hits (salp) > 1) 
 
1443
     {
 
1444
        newsalp = SortBlastHits (salp, msp->blast_threshold, msp->choice_blastfilter);
 
1445
     }
 
1446
     else newsalp = salp;
 
1447
 
 
1448
     if (number_hits (newsalp) > 1) 
 
1449
     {
 
1450
        if (newsalp != NULL && (msp->splicing)) 
 
1451
        {
 
1452
           headnode = CreateGraph (newsalp, msp->choice_blastfilter);
 
1453
           headnode = SimplifyGraph (headnode);
 
1454
           for(total=1, node=headnode; node!=NULL; node=node->next) 
 
1455
              total++; 
 
1456
           tab = (Int2Ptr)MemNew ((size_t)((total+2)*sizeof(Int2)));
 
1457
           for(j=0; j<total; j++)
 
1458
              tab[j] = -1;
 
1459
           vnp = NULL;
 
1460
           for (node=headnode; node!=NULL; node=node->next) {
 
1461
              vnp = traverseGraph (node, node, total, 0, tab, vnp, &change);
 
1462
           }
 
1463
           bestsol = get_solutions (vnp, headnode, total, SeqLocLen(slp2));
 
1464
           newsalp = link_solution (bestsol, headnode, total);
 
1465
        }
 
1466
        if (newsalp != NULL) 
 
1467
        {
 
1468
           if ((strand=SeqAlignStrand(newsalp, 0)) == Seq_strand_minus)
 
1469
           {
 
1470
              newsalp = SeqAlignListReverseStrand (newsalp);
 
1471
              newsalp = seqalign_reverse_sorting (newsalp);
 
1472
           }
 
1473
           newsalp = align_btwhits (newsalp, SeqLocId(slp1), SeqLocId(slp2), msp);
 
1474
        }
 
1475
     }
 
1476
     if (newsalp != NULL && msp->align_ends) {
 
1477
        newsalp = align_extrem (newsalp, slp1, slp2, msp);
 
1478
     }
 
1479
     if (delete_mash)
 
1480
        MemFree (msp);
 
1481
  }
 
1482
  return newsalp;
 
1483
}
 
1484
 
 
1485
 
 
1486
static SeqAlignPtr AlignAnyway (SeqLocPtr slp1, SeqLocPtr slp2, Boolean is_prot, MashPtr msp, Boolean message)
 
1487
{
 
1488
  SeqAlignPtr salp,
 
1489
              salpblast;
 
1490
  ValNodePtr  vnp;
 
1491
  Char        st1[50], st2[50];
 
1492
 
 
1493
  SeqIdWrite (SeqLocId(slp1), st1, PRINTID_FASTA_SHORT, 50);
 
1494
  SeqIdWrite (SeqLocId(slp2), st2, PRINTID_FASTA_SHORT, 50);
 
1495
  salpblast = BlastTwoSeqLocs (slp1, slp2, is_prot, msp); 
 
1496
  if (salpblast!=NULL) {
 
1497
     salp = BlastBandAlignTwoSeqLocs (slp1, slp2, salpblast, is_prot, msp);
 
1498
     if (salp!=NULL) {
 
1499
        if (message)
 
1500
           Message (MSG_OK, "%s - %s : Global alignment based on BLAST local similarity", st1, st2);
 
1501
        return salp;
 
1502
     }
 
1503
     SeqAlignSetFree (salpblast);
 
1504
  }
 
1505
  salp = BandAlignTwoSeqLocs (slp1, slp2, msp);
 
1506
  if (salp!=NULL) {
 
1507
     if (message)
 
1508
        Message (MSG_OK, "%s - %s : Global alignment using dynamic programing algorithm", st1, st2);
 
1509
     return salp;
 
1510
  }
 
1511
  vnp = NULL;
 
1512
  ValNodeAddPointer (&vnp, 0, slp1);
 
1513
  ValNodeAddPointer (&vnp, 0, slp2);
 
1514
  salp = SeqLocToFastaSeqAlign (vnp); 
 
1515
  ValNodeFreeType (&vnp, TypeEmpty);
 
1516
  if (salp!=NULL) {
 
1517
     if (message)
 
1518
        Message (MSG_OK, "Import sequence %s without alignment", st2);
 
1519
     return salp;
 
1520
  }
 
1521
  if (message)
 
1522
     ErrPostEx (SEV_WARNING, 0, 0, "No alignment");
 
1523
  return NULL;
 
1524
}
 
1525
 
 
1526
 
 
1527
static SeqAlignPtr AlignmRNA2genomicToSeqAlign (SeqLocPtr slp1, SeqLocPtr slp2, SeqAlignPtr salpblast, MashPtr msp)
 
1528
{
 
1529
  SeqAlignPtr salp=NULL;
 
1530
  Boolean     delete_mash = FALSE;
 
1531
  
 
1532
  if (salpblast != NULL) 
 
1533
  {
 
1534
     if (msp==NULL) {
 
1535
        msp = MashNew (FALSE);
 
1536
        delete_mash = TRUE;
 
1537
     }
 
1538
     if (msp!=NULL) {
 
1539
        msp->splicing = TRUE;
 
1540
        msp->choice_blastfilter = 2;
 
1541
        salp = BlastBandAlignTwoSeqLocs (slp1, slp2, salpblast, FALSE, msp);
 
1542
        if (delete_mash)
 
1543
           MemFree (msp);
 
1544
     }
 
1545
  }
 
1546
  return salp;
 
1547
}
 
1548
 
 
1549
/********************************************* 
 
1550
*** SeqLocListToSeqAlign 
 
1551
***    aligns the sequences from the SeqLocs list (sqloc_list) 
 
1552
***    returns a SeqAlign 
 
1553
***    Alignment methods:
 
1554
***      1: FASTA (assume that input is FASTA aligned )
 
1555
***      5: BandAlign (GlobalAlignTwoSeqLocs)
 
1556
**********************************************/
 
1557
 
 
1558
NLM_EXTERN SeqAlignPtr SeqLocListToSeqAlign (ValNodePtr sqloc_list, Int2 choice, Pointer param)
 
1559
{
 
1560
  ValNodePtr       vnp;
 
1561
  SeqLocPtr        master_slp, 
 
1562
                   slp=NULL;
 
1563
  SeqLocPtr        slptmp1,
 
1564
                   slptmp2;
 
1565
  SeqAlignPtr      salphead= NULL,
 
1566
                   salphead2 = NULL,
 
1567
                   salptmp = NULL,
 
1568
                   salpnew = NULL,
 
1569
                   salpblast = NULL,
 
1570
                   salpna = NULL,
 
1571
                   salpna2 = NULL;
 
1572
  SeqAnnotPtr      sapna;
 
1573
  MashPtr          msp;
 
1574
  ValNodePtr       framep = NULL;
 
1575
  Int2             seq_number = 0;
 
1576
  Uint1            frame;
 
1577
  Boolean          delete_mash = FALSE;
 
1578
  Boolean          is_prot = FALSE;
 
1579
 
 
1580
  BioseqPtr        bsp;
 
1581
 
 
1582
  if (sqloc_list == NULL) {
 
1583
      ErrPostEx(SEV_WARNING, 0, 0, "NULL SeqLocList");
 
1584
     return NULL;
 
1585
  }
 
1586
  if (choice==2 || choice==3 || choice==4 || choice==5)
 
1587
     return NULL;
 
1588
  if (choice == 1) 
 
1589
     return (SeqAlignPtr)SeqLocToFastaSeqAlign (sqloc_list);
 
1590
 
 
1591
  if((master_slp = (SeqLocPtr) sqloc_list->data.ptrvalue)==NULL) {
 
1592
     ErrPostEx(SEV_WARNING, 0, 0, "Can not read master sequence");
 
1593
     return NULL;
 
1594
  }
 
1595
  seq_number = 1;
 
1596
  msp = (MashPtr)param;
 
1597
  if (msp == NULL) {
 
1598
     bsp = BioseqLockById (SeqLocId(master_slp));
 
1599
     if (bsp) {
 
1600
        is_prot = ISA_aa(bsp->mol);
 
1601
        BioseqUnlock (bsp);
 
1602
     }
 
1603
     else
 
1604
        return NULL;
 
1605
     msp = MashNew (is_prot);
 
1606
     delete_mash = TRUE; 
 
1607
  }
 
1608
  if (msp == NULL)
 
1609
     return NULL;
 
1610
 
 
1611
  is_prot = (Boolean)(msp->is_prot || msp->translate_prot);
 
1612
  if (msp->translate_prot) {
 
1613
     slptmp1 = master_slp;
 
1614
     master_slp = TranslateSeqLoc(slptmp1, (Int2) Seq_code_ncbistdaa, &frame);
 
1615
     ValNodeAddInt (&framep, 1, (Int4)frame);
 
1616
  }
 
1617
/*   
 
1618
  if(msp->is_prot && msp->matrixname==NULL) msp->matrixname="BLOSUM62"; 
 
1619
*/
 
1620
  for (vnp = sqloc_list->next; vnp!=NULL; vnp = vnp->next)
 
1621
  {
 
1622
     salpnew = NULL;
 
1623
     slp = (SeqLocPtr) vnp->data.ptrvalue;
 
1624
     if ( slp != NULL ) 
 
1625
     {
 
1626
        if (msp->translate_prot) {
 
1627
           slptmp2 = slp;
 
1628
           slp = TranslateSeqLoc(slptmp2, (Int2) Seq_code_ncbistdaa, &frame);
 
1629
           ValNodeAddInt (&framep, 1, (Int4)frame);
 
1630
        }
 
1631
        if (choice == 2) {
 
1632
           ObjMgrSetHold ();
 
1633
/*
 
1634
COMMENT    salpnew = (SeqAlignPtr) SIM2ALN (master_slp, slp, 1); */
 
1635
           ObjMgrClearHold (); 
 
1636
        } 
 
1637
        else if (choice == 3) {
 
1638
           ObjMgrSetHold ();
 
1639
/*
 
1640
COMMENT    salpnew=(SeqAlignPtr)SIM3ALN_choice (master_slp, slp, (Int4) 100, TRUE, TRUE); */
 
1641
           ObjMgrClearHold (); 
 
1642
        } 
 
1643
        else if (choice == 4) {
 
1644
/*
 
1645
COMMENT    salpnew = (SeqAlignPtr) CSIM (master_slp, slp, 1, 0.00, NULL); */
 
1646
        } 
 
1647
        else if (choice == 5) {
 
1648
/*
 
1649
COMMENT    salpnew = SIM4ALN_choice (master_slp, slp, 1000, 8); */
 
1650
        } 
 
1651
        else if (choice == 6) {
 
1652
           salpblast = BlastTwoSeqLocs (master_slp, slp, is_prot, msp);
 
1653
           if (salpblast!=NULL) {
 
1654
              salpnew = BandAlignTwoSeqLocs (master_slp, slp, msp);
 
1655
           }
 
1656
           SeqAlignSetFree (salpblast);
 
1657
        } 
 
1658
        else if (choice == 7) {
 
1659
           salpnew = BlastTwoSeqLocs (master_slp, slp, is_prot, msp); 
 
1660
           msp->choice_blastfilter = 0;
 
1661
           salpnew = SortBlastHits (salpnew, msp->blast_threshold, msp->choice_blastfilter);
 
1662
        } 
 
1663
        else if (choice == 10) { 
 
1664
           salpblast = BlastTwoSeqLocs (master_slp, slp, is_prot, msp); 
 
1665
           if (salpblast!=NULL) 
 
1666
           {
 
1667
              salpnew = BlastBandAlignTwoSeqLocs (master_slp, slp, salpblast, is_prot, msp);
 
1668
              if (salpnew==NULL)
 
1669
                 SeqAlignSetFree (salpblast);
 
1670
           }
 
1671
        } 
 
1672
        else if (choice == 9) {
 
1673
           salpblast = BlastTwoSeqLocs (master_slp, slp, is_prot, msp); 
 
1674
           if (salpblast!=NULL) {
 
1675
              salpnew = AlignmRNA2genomicToSeqAlign (master_slp, slp, salpblast, msp);
 
1676
              if (salpnew==NULL)
 
1677
                 SeqAlignSetFree (salpblast);
 
1678
           }
 
1679
        } 
 
1680
        else if (choice == 8) {
 
1681
           salpblast = BlastTwoSeqLocs (master_slp, slp, is_prot, msp); 
 
1682
           if (salpblast!=NULL) {
 
1683
              salpnew = AlignmRNA2genomicToSeqAlign (master_slp, slp, salpblast, msp);
 
1684
              if (salpnew == NULL) 
 
1685
              {
 
1686
                 salpnew = BlastBandAlignTwoSeqLocs (master_slp, slp, salpblast, is_prot, msp);
 
1687
              }
 
1688
              if (salpnew==NULL)
 
1689
                 SeqAlignSetFree (salpblast);
 
1690
           }
 
1691
        }
 
1692
        else if (choice == 0)
 
1693
        {
 
1694
           if (seq_number<3)
 
1695
              salpnew = AlignAnyway (master_slp, slp, is_prot, msp, TRUE);
 
1696
           else
 
1697
              salpnew = AlignAnyway (master_slp, slp, is_prot, msp, FALSE);
 
1698
        }
 
1699
        if (salpnew != NULL) 
 
1700
        {
 
1701
           if (msp->translate_prot) {
 
1702
              master_slp = slptmp1;
 
1703
              sapna =SeqAnnotForSeqAlign (salpna);
 
1704
              sapna = aaSeqAnnot_to_dnaSeqAnnotFunc (&sapna, salpnew, sqloc_list, framep);
 
1705
              salpna = (SeqAlignPtr)sapna->data;
 
1706
              sapna->data=NULL;
 
1707
              SeqAnnotFree(sapna);
 
1708
           }
 
1709
           if (salphead == NULL) {
 
1710
              salphead = salpnew;
 
1711
              salptmp = salpnew;
 
1712
           }
 
1713
           else {
 
1714
              salptmp->next = salpnew;
 
1715
              salptmp = salptmp->next;
 
1716
           }
 
1717
           seq_number++;
 
1718
        }
 
1719
     } 
 
1720
     else {
 
1721
        ErrPostEx(SEV_WARNING, 0, 0, "Last SeqLoc was NULL");
 
1722
     }
 
1723
     
 
1724
  }
 
1725
  if (salphead != NULL) {
 
1726
     if (msp->map_align)
 
1727
        salphead = SeqAlignMapOnFirstSeq (salphead);
 
1728
     if (seq_number > 2 && msp->multidim) {
 
1729
        salphead2 = PairSeqAlign2MultiSeqAlign (salphead);
 
1730
        if (salphead2!=NULL)
 
1731
        {
 
1732
           salphead = SeqAlignSetFree (salphead);
 
1733
           salphead = salphead2;
 
1734
        }
 
1735
     }
 
1736
  }
 
1737
  if (salpna != NULL) {
 
1738
     if (seq_number > 2 && msp->multidim) {
 
1739
        salpna2 = PairSeqAlign2MultiSeqAlign (salpna);
 
1740
        if (salpna2!=NULL)
 
1741
        {
 
1742
           salpna = SeqAlignSetFree (salpna);
 
1743
           salpna = salpna2;
 
1744
        }
 
1745
     }
 
1746
     if (salpna !=NULL)
 
1747
        msp->transalp = salpna;
 
1748
  }
 
1749
  if (delete_mash)
 
1750
     MemFree (msp);
 
1751
  return salphead;
 
1752
}
 
1753
 
 
1754
 
 
1755
NLM_EXTERN SeqLocPtr AlignmRNA2genomic (BioseqPtr bsp1, BioseqPtr bsp2)
 
1756
{
 
1757
  ValNodePtr  vnp = NULL;
 
1758
  SeqLocPtr   slp = NULL,
 
1759
              slp1 = NULL,
 
1760
              slp2 = NULL;
 
1761
  SeqIdPtr    sip;
 
1762
  SeqAlignPtr salp,
 
1763
              salpblast;
 
1764
  MashPtr     msp = NULL;
 
1765
 
 
1766
  if (bsp1==NULL || bsp2==NULL)
 
1767
     return NULL;
 
1768
  if ((Boolean)ISA_aa(bsp1->mol) || (Boolean)ISA_aa(bsp2->mol))
 
1769
     return NULL;
 
1770
  sip = SeqIdFindBest (bsp1->id, 0);
 
1771
  if (sip==NULL)
 
1772
     return NULL;
 
1773
  slp1 = SeqLocIntNew (0, bsp1->length - 1, Seq_strand_plus, sip);
 
1774
  ValNodeAddPointer(&vnp, 0, (Pointer)slp1);
 
1775
  sip = SeqIdFindBest (bsp2->id, 0);
 
1776
  if (sip==NULL) 
 
1777
     return NULL;
 
1778
  slp2 = SeqLocIntNew (0, bsp2->length - 1, Seq_strand_plus, sip);
 
1779
  ValNodeAddPointer(&vnp, 0, (Pointer)slp2);
 
1780
  salpblast = BlastTwoSeqLocs (slp1, slp2, FALSE, msp);
 
1781
  if (salpblast!=NULL)
 
1782
  {
 
1783
     salp = AlignmRNA2genomicToSeqAlign (slp1, slp2, salpblast, msp);
 
1784
     SeqAlignSetFree (salpblast);
 
1785
     slp = SeqLocMixFromSeqAlign (salp, NULL);
 
1786
  }
 
1787
  return slp;
 
1788
}
 
1789
 
 
1790
NLM_EXTERN SeqAnnotPtr BlastBandAlignFromBlastSeqAlign (SeqAlignPtr salpblast, Boolean align_ends)
 
1791
{
 
1792
  BioseqPtr   bsp;
 
1793
  SeqAnnotPtr sap, sap2;
 
1794
  SeqAlignPtr salptmp,
 
1795
              salpnew,
 
1796
              salp_next,
 
1797
              salp_head=NULL,
 
1798
              salp_cp;
 
1799
  ValNodePtr  vnp_sameid=NULL,
 
1800
              vnp;
 
1801
  SeqLocPtr   master_slp, 
 
1802
              slp;
 
1803
  SeqIdPtr    sip,
 
1804
              siptmp;
 
1805
  Boolean     found=FALSE,
 
1806
              is_prot=FALSE;
 
1807
  MashPtr     msp = NULL;
 
1808
  Uint1       val;
 
1809
 
 
1810
  if (!salpblast)
 
1811
     return NULL;
 
1812
 
 
1813
  sap = SeqAnnotForSeqAlign (salpblast);
 
1814
 
 
1815
  sap2 = (SeqAnnotPtr) AsnIoMemCopy((Pointer)sap, 
 
1816
                (AsnReadFunc)SeqAnnotAsnRead, (AsnWriteFunc)SeqAnnotAsnWrite);
 
1817
 
 
1818
  salp_cp = (SeqAlignPtr)sap2->data;
 
1819
  sap->data=NULL;
 
1820
  SeqAnnotFree (sap);
 
1821
  sap2->data=NULL;
 
1822
  SeqAnnotFree (sap2);
 
1823
 
 
1824
  sip = SeqIdPtrFromSeqAlign (salp_cp);
 
1825
  bsp = BioseqLockById (sip);
 
1826
  if (bsp) {
 
1827
     is_prot = ISA_aa(bsp->mol);
 
1828
     master_slp = SeqLocIntNew (0, bsp->length-1, Seq_strand_plus, sip);
 
1829
     BioseqUnlock (bsp);
 
1830
  }
 
1831
  else return NULL;
 
1832
 
 
1833
  ValNodeAddPointer (&vnp_sameid, 0, (Pointer)salp_cp);
 
1834
  salp_next = salp_cp->next;
 
1835
  salp_cp->next = NULL;
 
1836
  salp_cp = salp_next;
 
1837
  while (salp_cp) {
 
1838
     sip = SeqIdPtrFromSeqAlign(salp_cp);
 
1839
     found = FALSE;
 
1840
     vnp = vnp_sameid;
 
1841
     while (vnp && !found) 
 
1842
     {
 
1843
        salptmp = (SeqAlignPtr)vnp->data.ptrvalue;
 
1844
        siptmp = SeqIdPtrFromSeqAlign (salptmp);
 
1845
        if ((val=SeqIdComp(sip->next, siptmp->next)) == SIC_YES)
 
1846
        {
 
1847
           while (salptmp->next!=NULL)
 
1848
              salptmp=salptmp->next;
 
1849
           salptmp->next = salp_cp;
 
1850
           salp_next = salp_cp->next;
 
1851
           salp_cp->next = NULL;
 
1852
           salp_cp = salp_next;
 
1853
           found = TRUE;
 
1854
        }
 
1855
        vnp=vnp->next;
 
1856
     }
 
1857
     if (!found) {
 
1858
        ValNodeAddPointer (&vnp_sameid, 0, (Pointer)salp_cp);
 
1859
        salp_next = salp_cp->next;
 
1860
        salp_cp->next = NULL;
 
1861
        salp_cp = salp_next;
 
1862
     }
 
1863
  }
 
1864
  if (vnp_sameid) {
 
1865
     msp = MashNew (is_prot);
 
1866
     msp->align_ends = align_ends;
 
1867
     for (vnp=vnp_sameid; vnp!=NULL; vnp=vnp->next) {
 
1868
        salptmp = (SeqAlignPtr)vnp->data.ptrvalue;
 
1869
        sip = SeqIdPtrFromSeqAlign (salptmp);
 
1870
        sip=sip->next;
 
1871
        bsp = BioseqLockById (sip);
 
1872
        if (bsp!= NULL)
 
1873
        {
 
1874
           slp = SeqLocIntNew (0, bsp->length-1, Seq_strand_plus, sip);
 
1875
           if (slp != NULL)
 
1876
           {
 
1877
              BioseqUnlock (bsp);
 
1878
              salpnew = BlastBandAlignTwoSeqLocs (master_slp, slp, salptmp, is_prot, msp);
 
1879
              if (salpnew) {
 
1880
                 if(salp_head == NULL)
 
1881
                    salp_head = salpnew;
 
1882
                 else {
 
1883
                    for(salptmp=salp_head; salptmp->next!=NULL; salptmp=salptmp->next)
 
1884
                       continue;
 
1885
                    salptmp->next=salpnew;
 
1886
                 }
 
1887
              }
 
1888
           }
 
1889
        }
 
1890
     }
 
1891
  }
 
1892
  if (salp_head==NULL)
 
1893
     return NULL;
 
1894
  return SeqAnnotForSeqAlign (salp_head);
 
1895
}
 
1896
 
 
1897
 
 
1898
/*******************************************************
 
1899
*** SeqEntryAlignToMasterFunc
 
1900
***    aligns the bioseqs that are present in a SeqEntry (sep)
 
1901
***    returns a SeqAlign
 
1902
*** SeqEntryToSeqAlign
 
1903
***    calls SeqEntryAlignToMaster
 
1904
***    returns a SeqAnnot
 
1905
******************************************************/
 
1906
 
 
1907
static SeqAlignPtr SeqEntryAlignToMasterFunc (SeqEntryPtr sep, SeqLocPtr master, Uint1 bsp_mol, Int2
 
1908
 method)
 
1909
{
 
1910
  ValNodePtr         vnp = NULL,
 
1911
                     vnp2 = NULL;
 
1912
  SeqAlignPtr        salp = NULL;
 
1913
  Int2               nb;
 
1914
 
 
1915
  vnp = SeqEntryToSeqLoc (sep, &nb, bsp_mol);
 
1916
  if (vnp != NULL) {
 
1917
     if (master != NULL)
 
1918
     {
 
1919
        ValNodeAddPointer (&vnp2, 0, master);
 
1920
        vnp2->next = vnp;
 
1921
        vnp = vnp2;
 
1922
     }
 
1923
     salp = SeqLocListToSeqAlign (vnp, method, NULL);
 
1924
  }
 
1925
  return salp;
 
1926
}
 
1927
 
 
1928
NLM_EXTERN SeqAnnotPtr LIBCALL SeqEntryToSeqAlign (SeqEntryPtr sep, Uint1 bsp_mol)
 
1929
{
 
1930
  SeqAlignPtr salp;
 
1931
  SeqAnnotPtr sap;
 
1932
 
 
1933
  salp = SeqEntryAlignToMasterFunc (sep, NULL, bsp_mol, PRG_ANYALIGN);
 
1934
  sap = SeqAnnotForSeqAlign (salp);
 
1935
/* this doesn't seem to produce a valid alignment and nobody knows what
 
1936
it's for anyway, so we'll comment it out and see who complains -- 8/8/01 SW
 
1937
  if (ISA_aa(bsp_mol)) {
 
1938
     sap = aaSeqAnnot_to_dnaSeqAnnotFunc (&sap, salp, NULL, NULL);
 
1939
  }
 
1940
*/
 
1941
  return sap;
 
1942
}