← Back to branch summary

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

« back to all changes in this revision

Viewing changes to tools/actutils.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:
tools/actutils.c
1
 
static char const rcsid[] = "$Id: actutils.c,v 6.46 2005/09/15 13:54:56 bollin Exp $";
 
1
static char const rcsid[] = "$Id: actutils.c,v 6.49 2006/01/04 18:52:34 bollin Exp $";
2
2
 
3
3
/* ===========================================================================
4
4
*
30
30
*
31
31
* Version Creation Date:   2/00
32
32
*
33
 
* $Revision: 6.46 $
 
33
* $Revision: 6.49 $
34
34
*
35
35
* File Description: utility functions for alignments
36
36
*
37
37
* Modifications:
38
38
* --------------------------------------------------------------------------
39
39
* $Log: actutils.c,v $
 
40
* Revision 6.49  2006/01/04 18:52:34  bollin
 
41
* changes to ACT_FindPiece to avoid memory leak
 
42
* changes to ACT_RemoveInconsistentAlnsFromSet to avoid global alignments that
 
43
* conflict by transposing elements
 
44
* changes to GetOldBlastAlignmentPiece to use a gap_open value that will not
 
45
* cause BLAST setup to fail
 
46
*
 
47
* Revision 6.48  2006/01/03 15:47:25  bollin
 
48
* allow alignment callback for ACT_FindPiece so that it can use the new BLAST
 
49
* library
 
50
*
 
51
* Revision 6.47  2005/12/19 14:32:54  bollin
 
52
* attempt to resolve conflicts between overlapping local alignments when
 
53
* constructing a global alignment - this resolved problems encountered when
 
54
* there are tandem repeats.
 
55
*
40
56
* Revision 6.46  2005/09/15 13:54:56  bollin
41
57
* don't create new entityID if bioseq already has one
42
58
*
189
205
#include <assert.h>
190
206
#include <viewmgr.h>
191
207
#include <alignmgr.h>
 
208
#include <salptool.h>
192
209
 
193
210
static void StateTableSearch (TextFsaPtr tbl, CharPtr txt, Int2Ptr state, Int4 pos, ACT_sitelistPtr PNTR asp_prev, ACT_sitelistPtr PNTR asp_head);
194
211
static Boolean am_isa_gap(Int4 start, Int4 prevstop, Uint1 strand);
1400
1417
   return 0;
1401
1418
}
1402
1419
 
 
1420
 
 
1421
/* spin structure: n1 = which alignment, n2 = start on first row, n3 =
 
1422
   alignment length on 1st row, n4 = start on 2nd row, n5 = 2nd strand */
 
1423
static void SetSpinValuesForSeqAlign (SeqAlignPtr salp, SQN_nPtr spin, Int4 list_pos, Int4 first_row)
 
1424
{
 
1425
   Int4             start;
 
1426
   Int4             stop;
 
1427
   Int4             strand;
 
1428
   
 
1429
   if (salp == NULL || spin == NULL || first_row > 2)
 
1430
   {
 
1431
      return;
 
1432
   }
 
1433
  
 
1434
   spin->n1 = list_pos;
 
1435
   AlnMgr2GetNthSeqRangeInSA(salp, first_row, &start, &stop);
 
1436
   spin->n3 = stop - start;
 
1437
   spin->n2 = start;
 
1438
   AlnMgr2GetNthSeqRangeInSA(salp, 3-first_row, &start, &stop);
 
1439
   spin->n4 = start;
 
1440
   strand = AlnMgr2GetNthStrand(salp, 3-first_row);
 
1441
   if (strand == Seq_strand_minus)
 
1442
      spin->n5 = -1;
 
1443
   else
 
1444
      spin->n5 = 1;
 
1445
  
 
1446
}
 
1447
 
 
1448
 
 
1449
static Boolean 
 
1450
ResolveSpinConflict 
 
1451
(SQN_nPtr    spin1,
 
1452
 SQN_nPtr    spin2, 
 
1453
 SeqAlignPtr salp, 
 
1454
 Int4        first_row)
 
1455
{
 
1456
   Boolean row1_conflict = FALSE, row2_conflict = FALSE;
 
1457
   Int4    left_aln_overlap = 0;
 
1458
   Int4    right_aln_overlap = 0;
 
1459
   Int4    spin1_row1_start, spin1_row1_stop;
 
1460
   Int4    spin2_row1_start, spin2_row1_stop;
 
1461
   Int4    spin1_row2_start, spin1_row2_stop;
 
1462
   Int4    spin2_row2_start, spin2_row2_stop;
 
1463
   Int4    aln_len;
 
1464
   Int4    strand;
 
1465
   Boolean resolved = FALSE;
 
1466
    
 
1467
   if (spin1 == NULL || spin2 == NULL || salp == NULL)
 
1468
   {
 
1469
      return FALSE;
 
1470
   }
 
1471
        
 
1472
   spin1_row1_start = spin1->n2;
 
1473
   spin1_row1_stop = spin1->n2 + spin1->n3 - 1;
 
1474
    
 
1475
   spin2_row1_start = spin2->n2;
 
1476
   spin2_row1_stop = spin2->n2 + spin2->n3 - 1;
 
1477
   
 
1478
   spin1_row2_start = spin1->n4;
 
1479
   spin1_row2_stop = spin1->n4 + spin1->n3 - 1;
 
1480
    
 
1481
   spin2_row2_start = spin2->n4;
 
1482
   spin2_row2_stop = spin2->n4 + spin2->n3 - 1;
 
1483
   
 
1484
   strand = spin1->n5;
 
1485
     
 
1486
   aln_len = SeqAlignLength (salp);
 
1487
      
 
1488
   /* NOTE - the spins were previously sorted, so it is impossible for
 
1489
    * spin2 to be longer than spin1
 
1490
    * Also, row 1 is required to be on the plus strand for both spins.
 
1491
    * We want to find the amount by which the alignment should be truncated,
 
1492
    * either on the left or on the right.  If the alignment needs to be truncated
 
1493
    * on both sides, there's probably something wrong.
 
1494
    * The alignment should not overlap on the first row or the second row.
 
1495
    */
 
1496
    
 
1497
   /* make sure we aren't aligning out of order */
 
1498
   if (strand == 1)
 
1499
   {
 
1500
     if ((spin1_row1_start > spin2_row1_start && spin1_row2_start < spin2_row2_start)
 
1501
         || (spin1_row1_start < spin2_row1_start && spin1_row2_start > spin2_row2_start))
 
1502
     {
 
1503
        return FALSE;
 
1504
     }
 
1505
   }
 
1506
   else
 
1507
   {
 
1508
     if ((spin1_row1_start > spin2_row1_start && spin1_row2_start > spin2_row2_start)
 
1509
         || (spin1_row1_start < spin2_row1_start && spin1_row2_start < spin2_row2_start))
 
1510
     {
 
1511
        return FALSE;
 
1512
     }
 
1513
   }
 
1514
    
 
1515
   /* check first for overlap on first row */
 
1516
   /* see if second is contained in first */
 
1517
   if (spin2_row1_start >= spin1_row1_start
 
1518
       && spin2_row1_stop <= spin1_row1_stop)
 
1519
   {
 
1520
      /*
 
1521
       * spin1 row 1: <---------->
 
1522
       * spin2 row 1:    <----->
 
1523
       */   
 
1524
      /* spin2 is contained in spin1 - no overlap could be removed */
 
1525
      return FALSE;
 
1526
   }
 
1527
   /* look for overlap on right */
 
1528
   if (spin2_row1_stop > spin1_row1_start
 
1529
       && spin2_row1_start < spin1_row1_start)
 
1530
   {
 
1531
      /*  
 
1532
       * spin1 row 1:    <----->
 
1533
       * spin2 row 1: <----->
 
1534
       */
 
1535
      right_aln_overlap = aln_len - AlnMgr2MapBioseqToSeqAlign (salp, spin1_row1_start, 1);
 
1536
   }
 
1537
   /* look for overlap on left */
 
1538
   if (spin2_row1_start < spin1_row1_stop
 
1539
       && spin2_row1_stop > spin1_row1_stop)
 
1540
   {
 
1541
      /*
 
1542
       * spin1 row 1: <----->
 
1543
       * spin2 row 1:    <----->
 
1544
       */   
 
1545
      left_aln_overlap = AlnMgr2MapBioseqToSeqAlign (salp, spin1_row1_stop, 1);
 
1546
   }
 
1547
   
 
1548
   /* check for overlap on second row */
 
1549
   /* see if second is contained in first */
 
1550
   if (spin2_row2_start >= spin1_row2_start
 
1551
       && spin2_row2_stop <= spin1_row2_stop)
 
1552
   {
 
1553
      /*
 
1554
       * spin1 row 2: <---------->
 
1555
       * spin2 row 2:    <----->
 
1556
       */   
 
1557
      /* spin2 is contained in spin1 - no overlap could be removed */
 
1558
      return FALSE;
 
1559
   }
 
1560
   
 
1561
   /* look for overlap on right */
 
1562
   if (spin2_row2_start < spin1_row2_start
 
1563
       && spin2_row2_stop > spin1_row2_start)
 
1564
   {
 
1565
      /*  
 
1566
       * spin1 row 2:    <----->
 
1567
       * spin2 row 2: <----->
 
1568
       */
 
1569
      if (strand == 1)
 
1570
      {
 
1571
        right_aln_overlap = MAX (right_aln_overlap, 
 
1572
                                 aln_len - AlnMgr2MapBioseqToSeqAlign (salp, spin1_row2_start, 1));
 
1573
      }
 
1574
      else
 
1575
      {
 
1576
        left_aln_overlap = MAX (left_aln_overlap, aln_len - AlnMgr2MapBioseqToSeqAlign (salp, spin1_row2_start, 1));
 
1577
      }
 
1578
   }
 
1579
   /* look for overlap on left */
 
1580
   if (spin2_row2_start < spin1_row2_stop
 
1581
       && spin2_row2_stop > spin1_row2_stop)
 
1582
   {
 
1583
      /*
 
1584
       * spin1 row 2: <----->
 
1585
       * spin2 row 2:    <----->
 
1586
       */   
 
1587
      if (strand == 1)
 
1588
      {
 
1589
         left_aln_overlap = MAX (left_aln_overlap, AlnMgr2MapBioseqToSeqAlign (salp, spin1_row2_stop, 1));
 
1590
      }
 
1591
      else
 
1592
      {
 
1593
        right_aln_overlap = MAX (right_aln_overlap, AlnMgr2MapBioseqToSeqAlign (salp, spin1_row2_stop, 1));
 
1594
      }
 
1595
   }   
 
1596
    
 
1597
   if (left_aln_overlap > 0 && right_aln_overlap == 0)
 
1598
   {
 
1599
      if (TruncateAlignment (salp, left_aln_overlap, TRUE))
 
1600
      {
 
1601
         SetSpinValuesForSeqAlign (salp, spin2, spin2->n1, first_row);
 
1602
         resolved = TRUE;
 
1603
      }
 
1604
   }
 
1605
   else if (left_aln_overlap == 0 && right_aln_overlap > 0)
 
1606
   {
 
1607
      if (TruncateAlignment (salp, right_aln_overlap, FALSE))
 
1608
      {
 
1609
         SetSpinValuesForSeqAlign (salp, spin2, spin2->n1, first_row);
 
1610
         resolved = TRUE;
 
1611
      }
 
1612
   }
 
1613
   
 
1614
   return resolved;   
 
1615
}
 
1616
 
 
1617
 
1403
1618
extern void ACT_RemoveInconsistentAlnsFromSet (SeqAlignPtr sap, Int4 fuzz, Int4 n)
1404
1619
{
1405
 
   AMAlignIndex2Ptr  amaip;
 
1620
   AMAlignIndex2Ptr amaip;
1406
1621
   Boolean          conflict;
1407
1622
   Int4             curr;
1408
1623
   Int4             i;
1411
1626
   SeqAlignPtr      salp_head;
1412
1627
   SeqAlignPtr      salp_prev;
1413
1628
   SQN_nPtr         PNTR spin;
1414
 
   Int4             start;
1415
 
   Int4             stop;
1416
1629
   Int4             strand;
1417
1630
 
1418
1631
   if (sap == NULL || sap->saip == NULL || sap->saip->indextype != INDEX_PARENT)
1441
1654
   for (i=0; i<amaip->numsaps; i++)
1442
1655
   {
1443
1656
      spin[i] = (SQN_nPtr)MemNew(sizeof(SQN_n));
1444
 
      salp = amaip->saps[i];
1445
 
      spin[i]->n1 = i;
1446
 
      AlnMgr2GetNthSeqRangeInSA(salp, n, &start, &stop);
1447
 
      spin[i]->n3 = stop - start;
1448
 
      spin[i]->n2 = start;
1449
 
      AlnMgr2GetNthSeqRangeInSA(salp, 3-n, &start, &stop);
1450
 
      spin[i]->n4 = start;
1451
 
      strand = AlnMgr2GetNthStrand(salp, 3-n);
1452
 
      if (strand == Seq_strand_minus)
1453
 
         spin[i]->n5 = -1;
1454
 
      else
1455
 
         spin[i]->n5 = 1;
 
1657
      SetSpinValuesForSeqAlign (amaip->saps[i], spin[i], i, n);
1456
1658
   }
1457
1659
   HeapSort((Pointer)spin, (size_t)(amaip->numsaps), sizeof(SQN_nPtr), ACT_CompareSpins);
1458
1660
   strand = spin[0]->n5;
1516
1718
                        conflict = TRUE;
1517
1719
                  }
1518
1720
               }
 
1721
               
 
1722
               /* make sure we aren't taking pieces out of order */
 
1723
               if (strand == 1)
 
1724
               {
 
1725
                  if ((spin[i]->n2 > spin[curr]->n2 && spin[i]->n4 < spin[curr]->n4)
 
1726
                      || (spin[i]->n2 < spin[curr]->n2 && spin[i]->n4 > spin[curr]->n4))
 
1727
                  {
 
1728
                    conflict = TRUE;
 
1729
                  }
 
1730
               }
 
1731
               else
 
1732
               {
 
1733
                  if ((spin[i]->n2 > spin[curr]->n2 && spin[i]->n4 > spin[curr]->n4)
 
1734
                      || (spin[i]->n2 < spin[curr]->n2 && spin[i]->n4 < spin[curr]->n4))
 
1735
                  {
 
1736
                    conflict = TRUE;
 
1737
                  }
 
1738
               }
 
1739
               
 
1740
               if (conflict && ResolveSpinConflict (spin[curr], spin[i], 
 
1741
                                                    amaip->saps[spin[i]->n1],
 
1742
                                                    n))
 
1743
               {
 
1744
                  conflict = FALSE;
 
1745
               }
 
1746
               
1519
1747
               if (conflict)
1520
1748
               {
1521
1749
                  salp = amaip->saps[spin[i]->n1];
1679
1907
   return sap;
1680
1908
}
1681
1909
 
1682
 
NLM_EXTERN SeqAlignPtr ACT_FindPiece(BioseqPtr bsp1, BioseqPtr bsp2, Int4 start1, Int4 stop1, Int4 start2, Int4 stop2, Uint1 strand, Int4 which_side)
 
1910
NLM_EXTERN SeqAlignPtr 
 
1911
ACT_FindPiece
 
1912
(BioseqPtr             bsp1,
 
1913
 BioseqPtr             bsp2, 
 
1914
 Int4                  start1, 
 
1915
 Int4                  stop1, 
 
1916
 Int4                  start2, 
 
1917
 Int4                  stop2, 
 
1918
 Uint1                 strand, 
 
1919
 Int4                  which_side,
 
1920
 GetAlignmentPieceFunc aln_piece_func)
1683
1921
{
1684
1922
   AMAlignIndex2Ptr      amaip;
1685
1923
   DenseSegPtr          dsp;
1688
1926
   Int4                 nstart2;
1689
1927
   Int4                 nstop1;
1690
1928
   Int4                 nstop2;
1691
 
   BLAST_OptionsBlkPtr  options;
1692
 
   CharPtr              program;
1693
1929
   SeqAlignPtr          sap;
1694
1930
   SeqAlignPtr          sap_head;
1695
1931
   SeqAlignPtr          sap_new;
1699
1935
 
1700
1936
   if (stop1 - start1 < 7 || stop2 - start2 < 7) /* can't do these by BLAST -- wordsize can't go that small */
1701
1937
      return NULL;
1702
 
   if (ISA_aa(bsp1->mol))
1703
 
   {
1704
 
      if (ISA_aa(bsp2->mol))
1705
 
         program = StringSave("blastp");
1706
 
      else
1707
 
         return NULL;
1708
 
   } else if (ISA_na(bsp1->mol))
1709
 
   {
1710
 
      if (ISA_na(bsp2->mol))
1711
 
         program = StringSave("blastn");
1712
 
      else
1713
 
         return NULL;
 
1938
   
 
1939
   if (aln_piece_func == NULL)
 
1940
   {
 
1941
      return NULL;
1714
1942
   }
1715
 
   options = BLASTOptionNew(program, TRUE);
1716
 
   options->gapped_calculation = TRUE;
1717
 
   options->expect_value = 10;
1718
 
   options->gap_x_dropoff_final = 100;
1719
 
   options->gap_open = 5;
1720
 
   options->gap_extend = 1;
1721
 
   options->penalty = -1;
1722
 
   options->wordsize = 7;
1723
1943
   slp1 = SeqLocIntNew(start1, stop1, Seq_strand_plus, bsp1->id);
1724
1944
   slp2 = SeqLocIntNew(start2, stop2, strand, bsp2->id);
1725
 
   sap = BlastTwoSequencesByLoc(slp1, slp2, program, options);
1726
 
   BLASTOptionDelete(options);
1727
 
   MemFree(program);
 
1945
   sap = aln_piece_func (slp1, slp2);
 
1946
 
 
1947
   SeqLocFree(slp1);
 
1948
   SeqLocFree(slp2);
 
1949
 
1728
1950
   if (sap == NULL)
1729
1951
   {
 
1952
      /* no real alignment, so create a "dummy" alignment that lines up at
 
1953
       * the requested side (left, right, or middle)
 
1954
       */
1730
1955
      sap = SeqAlignNew();
1731
1956
      dsp = DenseSegNew();
1732
1957
      dsp->numseg = 1;
1749
1974
      }
1750
1975
      sap->segs = (Pointer)dsp;
1751
1976
      sap->segtype = SAS_DENSEG;
1752
 
      return sap;
1753
 
   }
1754
 
   SeqLocFree(slp1);
1755
 
   SeqLocFree(slp2);
1756
 
   if (sap == NULL)
1757
 
      return NULL;
1758
 
   AlnMgr2IndexLite(sap);
1759
 
   ACT_RemoveInconsistentAlnsFromSet(sap, 20, 1);
1760
 
   amaip = (AMAlignIndex2Ptr)(sap->saip);
1761
 
   AlnMgr2SortAlnSetByNthRowPos(sap, 1);
1762
 
   ACT_GetNthSeqRangeInSASet(sap, 1, &nstart1, &nstop1);
1763
 
   ACT_GetNthSeqRangeInSASet(sap, 2, &nstart2, &nstop2);
1764
 
   strand = AlnMgr2GetNthStrand(amaip->saps[0], 2);
1765
 
   sap_head = NULL;
1766
 
   if (strand != Seq_strand_minus)
1767
 
   {
1768
 
      if (nstart1 > start1+20 && nstart2 > start2+20)
1769
 
      {
1770
 
         slp1 = SeqLocIntNew(start1, nstart1, Seq_strand_plus, bsp1->id);
1771
 
         slp2 = SeqLocIntNew(start2, nstart2, strand, bsp2->id);
1772
 
         sap_head = ACT_FindBestAlnByDotPlot(slp1, slp2);
1773
 
         SeqLocFree(slp1);
1774
 
         SeqLocFree(slp2);
1775
 
      }
1776
 
   } else
1777
 
   {
1778
 
      if (nstart1 > start1+20 && nstop2 < stop2 - 20)
1779
 
      {
1780
 
         slp1 = SeqLocIntNew(start1, nstart1, Seq_strand_plus, bsp1->id);
1781
 
         slp2 = SeqLocIntNew(nstop2, stop2, strand, bsp2->id);
1782
 
         sap_head = ACT_FindBestAlnByDotPlot(slp1, slp2);
1783
 
         SeqLocFree(slp1);
1784
 
         SeqLocFree(slp2);
1785
 
      }
1786
 
   }
1787
 
   for (i=0; i<amaip->numsaps-1; i++)
1788
 
   {
1789
 
      AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 1, NULL, &nstart1);
1790
 
      AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 1, &nstop1, NULL);
1791
 
      if (strand != Seq_strand_minus)
1792
 
      {
1793
 
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 2, NULL, &nstart2);
1794
 
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 2, &nstop2, NULL);
1795
 
      } else
1796
 
      {
1797
 
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 2, &nstop2, NULL);
1798
 
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 2, NULL, &nstart2);
1799
 
      }
1800
 
   }
1801
 
   ACT_GetNthSeqRangeInSASet(sap, 1, &nstart1, &nstop1);
1802
 
   ACT_GetNthSeqRangeInSASet(sap, 2, &nstart2, &nstop2);
1803
 
   sap_prev = sap_head;
1804
 
   if (sap_prev != NULL)
1805
 
   {
1806
 
      while (sap_prev->next != NULL)
1807
 
      {
1808
 
         sap_prev = sap_prev->next;
1809
 
      }
1810
 
   }
1811
 
   if (strand != Seq_strand_minus)
1812
 
   {
1813
 
      if (nstop1 < stop1-20 && nstop2 < stop2-20)  /* missing piece at the end */
1814
 
      {
1815
 
         slp1 = SeqLocIntNew(nstop1, stop1, Seq_strand_plus, bsp1->id);
1816
 
         slp2 = SeqLocIntNew(nstop2, stop2, strand, bsp2->id);
1817
 
         sap_new = ACT_FindBestAlnByDotPlot(slp1, slp2);
1818
 
         SeqLocFree(slp1);
1819
 
         SeqLocFree(slp2);
1820
 
         if (sap_new != NULL)
1821
 
         {
1822
 
            if (sap_head != NULL)
1823
 
            {
1824
 
               sap_prev->next = sap_new;
1825
 
               sap_prev = sap_new;
1826
 
            } else
1827
 
              sap_head = sap_prev = sap_new;
1828
 
         }
1829
 
      }
1830
 
   } else
1831
 
   {
1832
 
      if (nstop1 < stop1-20 && nstart2 > start2 + 20)
1833
 
      {
1834
 
         slp1 = SeqLocIntNew(nstop1, stop1, Seq_strand_plus, bsp1->id);
1835
 
         slp2 = SeqLocIntNew(start2, nstart2, strand, bsp2->id);
1836
 
         sap_new = ACT_FindBestAlnByDotPlot(slp1, slp2);
1837
 
         SeqLocFree(slp1);
1838
 
         SeqLocFree(slp2);
1839
 
         if (sap_new != NULL)
1840
 
         {
1841
 
            if (sap_head != NULL)
1842
 
            {
1843
 
               sap_prev->next = sap_new;
1844
 
               sap_prev = sap_new;
1845
 
            } else
1846
 
               sap_head = sap_prev = sap_new;
1847
 
         }
1848
 
      }
1849
 
   }
1850
 
   sap_new = (SeqAlignPtr)(sap->segs);
1851
 
   while (sap_new->next != NULL)
1852
 
   {
1853
 
      sap_new = sap_new->next;
1854
 
   }
1855
 
   sap_new->next = sap_head;
1856
 
   sap_head = (SeqAlignPtr)(sap->segs);
1857
 
   sap->segs = NULL;
1858
 
   SeqAlignFree(sap);
1859
 
   return sap_head;
 
1977
   }
 
1978
   else if (sap->next != NULL)
 
1979
   {
 
1980
      /* combine the alignments returned (if there are more than one) into
 
1981
       * a mini global alignment to cover the requested interval
 
1982
       */
 
1983
      AlnMgr2IndexLite(sap);
 
1984
      ACT_RemoveInconsistentAlnsFromSet(sap, 20, 1);
 
1985
      amaip = (AMAlignIndex2Ptr)(sap->saip);
 
1986
      AlnMgr2SortAlnSetByNthRowPos(sap, 1);
 
1987
      ACT_GetNthSeqRangeInSASet(sap, 1, &nstart1, &nstop1);
 
1988
      ACT_GetNthSeqRangeInSASet(sap, 2, &nstart2, &nstop2);
 
1989
      strand = AlnMgr2GetNthStrand(amaip->saps[0], 2);
 
1990
      sap_head = NULL;
 
1991
      if (strand != Seq_strand_minus)
 
1992
      {
 
1993
         if (nstart1 > start1+20 && nstart2 > start2+20)
 
1994
         {
 
1995
            slp1 = SeqLocIntNew(start1, nstart1, Seq_strand_plus, bsp1->id);
 
1996
            slp2 = SeqLocIntNew(start2, nstart2, strand, bsp2->id);
 
1997
            sap_head = ACT_FindBestAlnByDotPlot(slp1, slp2);
 
1998
            SeqLocFree(slp1);
 
1999
            SeqLocFree(slp2);
 
2000
         }
 
2001
      } 
 
2002
      else
 
2003
      {
 
2004
         if (nstart1 > start1+20 && nstop2 < stop2 - 20)
 
2005
         {
 
2006
            slp1 = SeqLocIntNew(start1, nstart1, Seq_strand_plus, bsp1->id);
 
2007
            slp2 = SeqLocIntNew(nstop2, stop2, strand, bsp2->id);
 
2008
            sap_head = ACT_FindBestAlnByDotPlot(slp1, slp2);
 
2009
            SeqLocFree(slp1);
 
2010
            SeqLocFree(slp2);
 
2011
         }
 
2012
      }
 
2013
      for (i=0; i<amaip->numsaps-1; i++)
 
2014
      {
 
2015
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 1, NULL, &nstart1);
 
2016
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 1, &nstop1, NULL);
 
2017
         if (strand != Seq_strand_minus)
 
2018
         {
 
2019
            AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 2, NULL, &nstart2);
 
2020
            AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 2, &nstop2, NULL);
 
2021
         } 
 
2022
         else
 
2023
         {
 
2024
            AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 2, &nstop2, NULL);
 
2025
            AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 2, NULL, &nstart2);
 
2026
         }
 
2027
      }
 
2028
      ACT_GetNthSeqRangeInSASet(sap, 1, &nstart1, &nstop1);
 
2029
      ACT_GetNthSeqRangeInSASet(sap, 2, &nstart2, &nstop2);
 
2030
      sap_prev = sap_head;
 
2031
      if (sap_prev != NULL)
 
2032
      {
 
2033
         while (sap_prev->next != NULL)
 
2034
         {
 
2035
            sap_prev = sap_prev->next;
 
2036
         }
 
2037
      }
 
2038
      if (strand != Seq_strand_minus)
 
2039
      {
 
2040
         if (nstop1 < stop1-20 && nstop2 < stop2-20)  /* missing piece at the end */
 
2041
         {
 
2042
            slp1 = SeqLocIntNew(nstop1, stop1, Seq_strand_plus, bsp1->id);
 
2043
            slp2 = SeqLocIntNew(nstop2, stop2, strand, bsp2->id);
 
2044
            sap_new = ACT_FindBestAlnByDotPlot(slp1, slp2);
 
2045
            SeqLocFree(slp1);
 
2046
            SeqLocFree(slp2);
 
2047
            if (sap_new != NULL)
 
2048
            {
 
2049
               if (sap_head != NULL)
 
2050
               {
 
2051
                  sap_prev->next = sap_new;
 
2052
                  sap_prev = sap_new;
 
2053
               } 
 
2054
               else
 
2055
                 sap_head = sap_prev = sap_new;
 
2056
            }
 
2057
         }
 
2058
      } 
 
2059
      else
 
2060
      {
 
2061
         if (nstop1 < stop1-20 && nstart2 > start2 + 20)
 
2062
         {
 
2063
            slp1 = SeqLocIntNew(nstop1, stop1, Seq_strand_plus, bsp1->id);
 
2064
            slp2 = SeqLocIntNew(start2, nstart2, strand, bsp2->id);
 
2065
            sap_new = ACT_FindBestAlnByDotPlot(slp1, slp2);
 
2066
            SeqLocFree(slp1);
 
2067
            SeqLocFree(slp2);
 
2068
            if (sap_new != NULL)
 
2069
            {
 
2070
               if (sap_head != NULL)
 
2071
               {
 
2072
                  sap_prev->next = sap_new;
 
2073
                  sap_prev = sap_new;
 
2074
               } 
 
2075
               else
 
2076
                  sap_head = sap_prev = sap_new;
 
2077
            }
 
2078
         }
 
2079
      }
 
2080
      sap_new = (SeqAlignPtr)(sap->segs);
 
2081
      while (sap_new->next != NULL)
 
2082
      {
 
2083
         sap_new = sap_new->next;
 
2084
      }
 
2085
      sap_new->next = sap_head;
 
2086
      sap_head = (SeqAlignPtr)(sap->segs);
 
2087
      sap->segs = NULL;
 
2088
      SeqAlignFree(sap);
 
2089
      sap = sap_head;
 
2090
   }
 
2091
   return sap;
1860
2092
}
1861
2093
 
1862
2094
static void ACT_ExtendAlnRight(SeqAlignPtr sap, Int4 which_row, Int4 start, Int4 stop)
2646
2878
   return sap;
2647
2879
}
2648
2880
 
2649
 
NLM_EXTERN SeqAlignPtr Sqn_GlobalAlign2SeqEx (BioseqPtr bsp1, BioseqPtr bsp2, BoolPtr revcomp, GetAlignmentFunc aln_func)
 
2881
static SeqAlignPtr LIBCALLBACK GetOldBlastAlignmentPiece (SeqLocPtr slp1, SeqLocPtr slp2)
 
2882
{
 
2883
   BioseqPtr            bsp1, bsp2;
 
2884
   BLAST_OptionsBlkPtr  options;
 
2885
   CharPtr              program;
 
2886
   SeqAlignPtr          sap = NULL;
 
2887
   
 
2888
   if (slp1 == NULL || slp2 == NULL)
 
2889
   { 
 
2890
      return NULL;
 
2891
   }
 
2892
   
 
2893
   bsp1 = BioseqFindFromSeqLoc (slp1);
 
2894
   bsp2 = BioseqFindFromSeqLoc (slp2);
 
2895
   
 
2896
   if (bsp1 == NULL || bsp2 == NULL)
 
2897
   {
 
2898
      return NULL;
 
2899
   }
 
2900
   
 
2901
   if (ISA_aa(bsp1->mol))
 
2902
   {
 
2903
      if (ISA_aa(bsp2->mol))
 
2904
         program = StringSave("blastp");
 
2905
      else
 
2906
         return NULL;
 
2907
   } else if (ISA_na(bsp1->mol))
 
2908
   {
 
2909
      if (ISA_na(bsp2->mol))
 
2910
         program = StringSave("blastn");
 
2911
      else
 
2912
         return NULL;
 
2913
   }
 
2914
   options = BLASTOptionNew(program, TRUE);
 
2915
   options->gapped_calculation = TRUE;
 
2916
   options->expect_value = 10;
 
2917
   options->gap_x_dropoff_final = 100;
 
2918
   options->gap_open = 4;
 
2919
   options->gap_extend = 1;
 
2920
   options->penalty = -1;
 
2921
   options->wordsize = 7;
 
2922
   sap = BlastTwoSequencesByLoc(slp1, slp2, program, options);
 
2923
   BLASTOptionDelete(options);
 
2924
   MemFree(program);
 
2925
   return sap;
 
2926
}
 
2927
 
 
2928
NLM_EXTERN SeqAlignPtr Sqn_GlobalAlign2SeqEx (BioseqPtr bsp1, BioseqPtr bsp2, BoolPtr revcomp, GetAlignmentFunc aln_func, GetAlignmentPieceFunc aln_piece_func)
2650
2929
{
2651
2930
   AMAlignIndex2Ptr     amaip;
2652
2931
   Int4                 i;
2709
2988
   sap_head = NULL;
2710
2989
 
2711
2990
   if (start1 > 6 && start1 < extnd)
2712
 
      sap_head = sap_prev = ACT_FindPiece(bsp1, bsp2, MAX(start1-SQN_WINDOW, 0), start1, MAX(start1-SQN_WINDOW, 0), start2, strand, SQN_LEFT);
 
2991
      sap_head = sap_prev = ACT_FindPiece(bsp1, bsp2, MAX(start1-SQN_WINDOW, 0), start1, MAX(start1-SQN_WINDOW, 0), start2, strand, SQN_LEFT, aln_piece_func);
2713
2992
   else if (start1 > 0 && start1 < extnd)
2714
2993
      SQN_ExtendAlnAlg(amaip->saps[0], start1, SQN_LEFT, Seq_strand_plus);
2715
2994
   ACT_GetNthSeqRangeInSASet(sap, 1, &start1, &stop1);
2716
2995
   ACT_GetNthSeqRangeInSASet(sap, 2, &start2, &stop2);
2717
2996
   if (start2 > 6 && start2 < extnd)
2718
2997
   {
2719
 
      sap_new = ACT_FindPiece(bsp2, bsp1, MAX(start2-SQN_WINDOW, 0), start2, MAX(start1-SQN_WINDOW, 0), start1, strand, SQN_LEFT);
 
2998
      sap_new = ACT_FindPiece(bsp2, bsp1, MAX(start2-SQN_WINDOW, 0), start2, MAX(start1-SQN_WINDOW, 0), start1, strand, SQN_LEFT, aln_piece_func);
2720
2999
      if (sap_new != NULL)
2721
3000
         SPI_flip_sa_list(sap_new);
2722
3001
      if (sap_head != NULL)
2740
3019
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i], 2, &stop2, NULL);
2741
3020
         AlnMgr2GetNthSeqRangeInSA(amaip->saps[i+1], 2, NULL, &start2);
2742
3021
      }
2743
 
      sap_new = ACT_FindPiece(bsp1, bsp2, start1, stop1, start2, stop2, strand, SQN_MIDDLE);
 
3022
      sap_new = ACT_FindPiece(bsp1, bsp2, start1, stop1, start2, stop2, strand, SQN_MIDDLE, aln_piece_func);
2744
3023
      if (sap_head)
2745
3024
      {
2746
3025
         sap_prev->next = sap_new;
2753
3032
   ACT_GetNthSeqRangeInSASet(sap, 2, &start2, &stop2);
2754
3033
   if (bsp1->length-stop1 > 6 && bsp1->length-stop1 < extnd)
2755
3034
   {
2756
 
      sap_new = ACT_FindPiece(bsp1, bsp2, stop1, MIN(bsp1->length-1, stop1 + SQN_WINDOW), stop2, MIN(bsp2->length-1, stop2+SQN_WINDOW), strand, SQN_RIGHT);
 
3035
      sap_new = ACT_FindPiece(bsp1, bsp2, stop1, MIN(bsp1->length-1, stop1 + SQN_WINDOW), stop2, MIN(bsp2->length-1, stop2+SQN_WINDOW), strand, SQN_RIGHT, aln_piece_func);
2757
3036
      if (sap_new != NULL)
2758
3037
      {
2759
3038
         if (sap_head != NULL)
2769
3048
   ACT_GetNthSeqRangeInSASet(sap, 2, &start2, &stop2);
2770
3049
   if (bsp2->length-stop2 > 6 && bsp2->length-stop2 < extnd)
2771
3050
   {
2772
 
      sap_new = ACT_FindPiece(bsp2, bsp1, stop2, MIN(bsp2->length-1, stop2 + SQN_WINDOW), stop1, MIN(bsp1->length-1, stop1+SQN_WINDOW), strand, SQN_RIGHT);
 
3051
      sap_new = ACT_FindPiece(bsp2, bsp1, stop2, MIN(bsp2->length-1, stop2 + SQN_WINDOW), stop1, MIN(bsp1->length-1, stop1+SQN_WINDOW), strand, SQN_RIGHT, aln_piece_func);
2773
3052
      if (sap_new != NULL)
2774
3053
      {
2775
3054
         SPI_flip_sa_list(sap_new);
2798
3077
 
2799
3078
NLM_EXTERN SeqAlignPtr Sqn_GlobalAlign2Seq (BioseqPtr bsp1, BioseqPtr bsp2, BoolPtr revcomp)
2800
3079
{
2801
 
  return Sqn_GlobalAlign2SeqEx (bsp1, bsp2, revcomp, GetOldBlastAlignment);
 
3080
  return Sqn_GlobalAlign2SeqEx (bsp1, bsp2, revcomp, GetOldBlastAlignment, GetOldBlastAlignmentPiece);
2802
3081
}
2803
3082
 
2804
3083