← Back to branch summary

~ubuntu-branches/ubuntu/precise/ncbi-tools6/precise

« back to all changes in this revision

Viewing changes to biostruc/cdd/cddutil.c

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
biostruc/cdd/cddutil.c
1
 
/* $Id: cddutil.c,v 1.45 2001/11/15 15:35:13 kans Exp $
 
1
/* $Id: cddutil.c,v 1.92 2004/06/22 14:20:19 camacho Exp $
2
2
*===========================================================================
3
3
*
4
4
*                            PUBLIC DOMAIN NOTICE
29
29
*
30
30
* Initial Version Creation Date: 10/18/1999
31
31
*
32
 
* $Revision: 1.45 $
 
32
* $Revision: 1.92 $
33
33
*
34
34
* File Description: CDD utility routines
35
35
*
36
36
* Modifications:
37
37
* --------------------------------------------------------------------------
38
38
* $Log: cddutil.c,v $
 
39
* Revision 1.92  2004/06/22 14:20:19  camacho
 
40
* Changed invocation of posFreqsToMatrix to conform with new signature.
 
41
*
 
42
* Revision 1.91  2003/12/09 22:21:35  bauer
 
43
* added CddCountResTypes
 
44
*
 
45
* Revision 1.90  2003/12/09 18:48:05  bauer
 
46
* commented out deprecated fields in BlastKarlinBlkCopy
 
47
*
 
48
* Revision 1.89  2003/10/07 21:19:39  bauer
 
49
* made CddMasterIs3D more general
 
50
*
 
51
* Revision 1.88  2003/08/25 19:09:47  bauer
 
52
* added SeqAlignReadFromFile
 
53
*
 
54
* Revision 1.87  2003/05/22 18:39:14  bauer
 
55
* list txids in Cdd XML dumper
 
56
*
 
57
* Revision 1.86  2003/05/21 17:25:21  bauer
 
58
* optional ObjMgr in CddReadDBGetBioseq
 
59
*
 
60
* Revision 1.85  2003/05/09 20:38:12  bauer
 
61
* report invalid intervals in CddRelocateSeqLoc
 
62
*
 
63
* Revision 1.84  2003/05/08 14:40:09  bauer
 
64
* bugfix in CddConsensus
 
65
*
 
66
* Revision 1.83  2003/04/25 17:17:32  thiessen
 
67
* fix return value type
 
68
*
 
69
* Revision 1.82  2003/04/25 14:36:20  bauer
 
70
* impalaScaling now returns boolean value
 
71
*
 
72
* Revision 1.81  2003/02/06 21:04:27  bauer
 
73
* fixed bug in reindexing to consensus
 
74
*
 
75
* Revision 1.80  2002/12/03 14:36:31  bauer
 
76
* added CddMSLMixedToMSLDenDiag
 
77
*
 
78
* Revision 1.79  2002/11/22 21:35:23  bauer
 
79
* added SeqAnnotReadFromFile and preservation of scores in DenseSeg to DenseDiag conversion
 
80
*
 
81
* Revision 1.78  2002/10/22 14:54:27  bauer
 
82
* fixed date format in XML dumper
 
83
*
 
84
* Revision 1.77  2002/10/10 20:38:19  bauer
 
85
* changes to accomodate new spec items
 
86
* - old-root node
 
87
* - curation-status
 
88
*
 
89
* Revision 1.76  2002/10/02 17:32:21  bauer
 
90
* avoid merging blocks when reindexing alignments
 
91
*
 
92
* Revision 1.75  2002/08/17 11:55:07  bauer
 
93
* backed out changes
 
94
*
 
95
* Revision 1.74  2002/08/16 19:51:45  bauer
 
96
* added Ben's CddSrvGetStyle2
 
97
*
 
98
* Revision 1.73  2002/08/06 12:54:26  bauer
 
99
* fixes to accomodate COGs
 
100
*
 
101
* Revision 1.72  2002/07/31 17:15:29  kans
 
102
* added prototype for BlastDefLineSetFree for Mac compiler
 
103
*
 
104
* Revision 1.71  2002/07/31 14:58:58  bauer
 
105
* BLAST DB Sequence Retrieval
 
106
*
 
107
* Revision 1.70  2002/07/17 18:54:10  bauer
 
108
* CddFeaturesAreConsistent now returns explicit error messages
 
109
*
 
110
* Revision 1.69  2002/07/11 14:43:43  bauer
 
111
* column 21(X) is now always -1 in PSSMs
 
112
*
 
113
* Revision 1.68  2002/07/10 22:51:10  bauer
 
114
* changed score for 26th pssm column to match whats done by makemat
 
115
*
 
116
* Revision 1.67  2002/07/10 22:06:07  bauer
 
117
* replaced posScaling by impalaScaling
 
118
*
 
119
* Revision 1.66  2002/07/10 15:34:31  bauer
 
120
* made SipIsConsensus public
 
121
*
 
122
* Revision 1.65  2002/07/09 21:12:39  bauer
 
123
* added CddDenDiagCposComp2KBP to return Karlin-Altschul parameters together with PSSM
 
124
*
 
125
* Revision 1.64  2002/07/05 21:09:25  bauer
 
126
* added GetAlignmentSize
 
127
*
 
128
* Revision 1.63  2002/06/12 15:22:51  bauer
 
129
* 6/11/02 update
 
130
*
 
131
* Revision 1.62  2002/05/31 17:54:31  thiessen
 
132
* fix for read-only string literal (e.g. on Mac)
 
133
*
 
134
* Revision 1.61  2002/05/24 17:49:01  bauer
 
135
* Unlock Bioseqs again in SeqAlignInform
 
136
*
 
137
* Revision 1.60  2002/05/16 22:37:49  bauer
 
138
* free seqalign in CddExpAlignToSeqAlign if empty
 
139
*
 
140
* Revision 1.59  2002/05/06 16:59:50  bauer
 
141
* remove blanks in inferred CD short names
 
142
*
 
143
* Revision 1.58  2002/04/25 14:30:19  bauer
 
144
* fixed CddFindMMDBIdInBioseq
 
145
*
 
146
* Revision 1.57  2002/04/22 16:37:31  bauer
 
147
* added check for missing structure alignments
 
148
*
 
149
* Revision 1.56  2002/04/18 21:00:26  bauer
 
150
* added check CddFeaturesAreConsistent
 
151
*
 
152
* Revision 1.55  2002/04/12 14:02:43  bauer
 
153
* added update_date case to CddAssignDescr
 
154
*
 
155
* Revision 1.54  2002/04/11 14:34:25  bauer
 
156
* added CddRemoveConsensus
 
157
*
 
158
* Revision 1.53  2002/03/28 15:55:00  bauer
 
159
* fixed bug in CddFindBioseqInSeqEntry, thanks to Ben
 
160
*
 
161
* Revision 1.52  2002/02/20 22:27:28  bauer
 
162
* utility functions for the CD-Server
 
163
*
 
164
* Revision 1.51  2002/02/12 23:00:46  bauer
 
165
* added missing break in CddRelocateSeqLoc
 
166
*
 
167
* Revision 1.50  2002/02/06 19:35:37  bauer
 
168
* some more CddGet.. functionality
 
169
*
 
170
* Revision 1.49  2002/02/05 23:15:40  bauer
 
171
* added a few CddGet.. utility functions, changes to CddAddDescr allow to extend scrapbook line by line
 
172
*
 
173
* Revision 1.48  2002/01/28 14:11:29  bauer
 
174
* add score to entrez neighbor lists
 
175
*
 
176
* Revision 1.47  2002/01/05 14:49:44  bauer
 
177
* made SeqAlignDup a local function
 
178
*
 
179
* Revision 1.46  2002/01/04 19:46:56  bauer
 
180
* added functions to interconvert PSSM-Ids and accessions
 
181
*
39
182
* Revision 1.45  2001/11/15 15:35:13  kans
40
183
* changed strdup to StringSave for Mac
41
184
*
188
331
#include <objmime.h>
189
332
#include <strimprt.h>
190
333
#include <cdd.h>
191
 
/* #include <bjcdd.h> */
 
334
/* #include <objcdd.h> */
 
335
#include <readdb.h>
 
336
#include <fdlobj.h>
192
337
#include <cddutil.h>
193
338
#include <edutil.h>
194
339
#include <posit.h>
231
376
  CddPtr   pcdd;
232
377
  
233
378
  if (bBin) {
234
 
     aip = AsnIoOpen(cFile,"rb");
235
 
     pcdd = CddAsnRead(aip,NULL);
 
379
    aip = AsnIoOpen(cFile,"rb");
236
380
  } else {
237
 
     aip = AsnIoOpen(cFile,"r");
238
 
     pcdd = CddAsnRead(aip,NULL);
 
381
    aip = AsnIoOpen(cFile,"r");
239
382
  }
240
 
 
 
383
  pcdd = CddAsnRead(aip,NULL);
241
384
  AsnIoClose(aip);
242
385
  return(pcdd);
243
386
}
244
387
 
 
388
SeqAnnotPtr LIBCALL SeqAnnotReadFromFile(CharPtr cFile, Boolean bBin)
 
389
{
 
390
  AsnIoPtr    aip = NULL;
 
391
  SeqAnnotPtr sap;
 
392
 
 
393
  if (bBin) {
 
394
    aip = AsnIoOpen(cFile,"rb");
 
395
  } else {
 
396
    aip = AsnIoOpen(cFile,"r");
 
397
  }
 
398
  sap = SeqAnnotAsnRead(aip, NULL);
 
399
  AsnIoClose(aip);
 
400
  return(sap);
 
401
}
 
402
 
 
403
SeqAlignPtr LIBCALL SeqAlignReadFromFile(CharPtr cFile, Boolean bBin)
 
404
{
 
405
  AsnIoPtr    aip = NULL;
 
406
  SeqAlignPtr salp;
 
407
 
 
408
  if (bBin) {
 
409
    aip = AsnIoOpen(cFile,"rb");
 
410
  } else {
 
411
    aip = AsnIoOpen(cFile,"r");
 
412
  }
 
413
  salp = SeqAlignAsnRead(aip, NULL);
 
414
  AsnIoClose(aip);
 
415
  return(salp);
 
416
}
 
417
 
245
418
/*---------------------------------------------------------------------------*/
246
419
/*---------------------------------------------------------------------------*/
247
420
/* Cdd-tree asn1 reader and writer wrappers                                  */
288
461
/*---------------------------------------------------------------------------*/
289
462
void LIBCALL CddAssignDescr(CddPtr pcdd, Pointer pThis, Int4 iWhat, Int4 iIval)
290
463
{
291
 
  ValNodePtr  vnp;
 
464
  ValNodePtr  vnp, vnp2;
292
465
 
293
466
  vnp = ValNodeNew(NULL);
294
467
  vnp->choice = iWhat;
306
479
      vnp->data.ptrvalue = (ValNodePtr) pThis;
307
480
      break;
308
481
    case CddDescr_create_date:
 
482
    case CddDescr_update_date:
309
483
      vnp->data.ptrvalue = (DatePtr) pThis;
310
484
      break;
311
485
    case CddDescr_tax_source:
317
491
    case CddDescr_status:
318
492
      vnp->data.intvalue = (Int4) iIval;
319
493
      break;
 
494
    case CddDescr_curation_status:
 
495
      vnp->data.intvalue = (Int4) iIval;
 
496
      break;
 
497
    case CddDescr_old_root:
 
498
      vnp->data.ptrvalue = (ValNodePtr) pThis;
 
499
      break;
320
500
    case CddDescr_scrapbook:
321
 
      vnp->data.ptrvalue = (ValNodePtr) pThis;
 
501
      vnp = ValNodeExtractList(&pcdd->description, CddDescr_scrapbook);
 
502
      vnp2 = ValNodeCopyStr(NULL,0,(CharPtr)pThis);
 
503
      if(vnp) {
 
504
        if(vnp->next) ErrPostEx(SEV_WARNING,0,0,"Multiple scrapbooks");
 
505
        ValNodeLink((ValNodePtr *)&vnp->data.ptrvalue, vnp2);
 
506
      } else ValNodeAddPointer(&vnp, CddDescr_scrapbook, vnp2);
 
507
      break;
322
508
    default:
323
509
      vnp->data.ptrvalue = pThis;
324
510
      break;
326
512
  ValNodeLink(&(pcdd->description),vnp);
327
513
}
328
514
 
 
515
 
 
516
/*---------------------------------------------------------------------------*/
 
517
/*---------------------------------------------------------------------------*/
 
518
/* remove an item from the CD description. Works for a limited set of types. */
 
519
/* status and repeats are removed no matter what. For comments, references   */
 
520
/* etc. the content of the description is compared to the input              */
 
521
/*---------------------------------------------------------------------------*/
 
522
/*---------------------------------------------------------------------------*/
 
523
Boolean LIBCALL CddKillDescr(CddPtr pcdd, Pointer pThis, Int4 iWhat, Int4 iIval) 
 
524
{
 
525
  ValNodePtr   vnp, vnpold = NULL, vnpthis = NULL, vnpbefore = NULL, pub;
 
526
  
 
527
  vnp = pcdd->description;
 
528
  while (vnp) {
 
529
    if (vnp->choice == iWhat) {
 
530
      switch (iWhat) {
 
531
        case CddDescr_comment:
 
532
        case CddDescr_othername:
 
533
        case CddDescr_category:
 
534
        case CddDescr_source:
 
535
          if (Nlm_StrCmp((CharPtr)pThis,(CharPtr)vnp->data.ptrvalue) == 0) {
 
536
            if (!vnpthis) {
 
537
              vnpthis = vnp;
 
538
              vnpbefore = vnpold;
 
539
            }
 
540
          }
 
541
          break;
 
542
        case CddDescr_status:
 
543
        case CddDescr_curation_status:
 
544
        case CddDescr_old_root:
 
545
        case CddDescr_repeats:
 
546
          if (!vnpthis) {
 
547
            vnpthis = vnp;
 
548
            vnpbefore = vnpold;
 
549
          }
 
550
          break;
 
551
        case CddDescr_reference:
 
552
          pub = vnp->data.ptrvalue;
 
553
          if ((pThis && pThis == pub) || 
 
554
              (iIval>0 && pub->choice==PUB_PMid && pub->data.intvalue==iIval)) {
 
555
            if (!vnpthis) {
 
556
              vnpthis = vnp;
 
557
              vnpbefore = vnpold;
 
558
            }
 
559
          }
 
560
        default:
 
561
          break;
 
562
      }    
 
563
    }
 
564
    vnpold = vnp;
 
565
    vnp = vnp->next;
 
566
  }
 
567
  if (vnpthis) {
 
568
    if (vnpbefore) {
 
569
      vnpbefore->next = vnpthis->next;    
 
570
    } else {
 
571
      pcdd->description = vnpthis->next;
 
572
    }
 
573
    return TRUE;
 
574
  } else return FALSE;
 
575
}
 
576
 
 
577
/*---------------------------------------------------------------------------*/
 
578
/*---------------------------------------------------------------------------*/
 
579
/* series of functions to extract info from CDs                              */
 
580
/*---------------------------------------------------------------------------*/
 
581
/*---------------------------------------------------------------------------*/
 
582
Int4 LIBCALL CddGetVersion(CddPtr pcdd) 
 
583
{
 
584
  CddIdPtr    cid;
 
585
  GlobalIdPtr pGid;
 
586
 
 
587
  cid = pcdd->id; while (cid) {
 
588
    if (cid->choice == CddId_gid) {
 
589
      pGid = cid->data.ptrvalue;
 
590
      return pGid->version;
 
591
    }
 
592
    cid = cid->next;
 
593
  }
 
594
  return 0;
 
595
}
 
596
 
 
597
OrgRefPtr LIBCALL CddGetOrgRef(CddPtr pcdd)
 
598
{
 
599
  CddDescrPtr   pcdsc;
 
600
 
 
601
  pcdsc = pcdd->description;
 
602
  while (pcdsc) {
 
603
    if (pcdsc->choice == CddDescr_tax_source) {
 
604
      return (OrgRefPtr) pcdsc->data.ptrvalue;
 
605
    }
 
606
    pcdsc = pcdsc->next;
 
607
  }
 
608
  return NULL;
 
609
}
 
610
 
 
611
Int4 LIBCALL CddGetPssmId(CddPtr pcdd) 
 
612
{
 
613
  CddIdPtr cid;
 
614
  
 
615
  cid = pcdd->id; while (cid) {
 
616
    if (cid->choice == CddId_uid) {
 
617
      return cid->data.intvalue;
 
618
    }
 
619
    cid = cid->next;
 
620
  }
 
621
  return 0;
 
622
}
 
623
 
 
624
Int4 LIBCALL CddGetPmIds(CddPtr pcdd, Int4Ptr iPMids)
 
625
{
 
626
  Int4          count = 0;
 
627
  CddDescrPtr   pcdsc;
 
628
  ValNodePtr    pub;
 
629
  Boolean       bRefOpen = FALSE;
 
630
 
 
631
  if (iPMids) iPMids = MemNew(250*sizeof(Int4));
 
632
  pcdsc = pcdd->description;
 
633
  while (pcdsc) {
 
634
    if (pcdsc->choice == CddDescr_reference) {
 
635
      pub = pcdsc->data.ptrvalue;
 
636
      if (pub->choice == PUB_PMid) {
 
637
        iPMids[count] = pub->data.intvalue;
 
638
        count++;
 
639
        if (count >= 250) return count;
 
640
      }
 
641
    }
 
642
    pcdsc = pcdsc->next;
 
643
  }
 
644
  return count;
 
645
}
 
646
 
 
647
 
 
648
CharPtr LIBCALL CddGetDescr(CddPtr pcdd)
 
649
{
 
650
  CddDescrPtr   pcdsc;
 
651
 
 
652
  pcdsc = pcdd->description;
 
653
  while (pcdsc) {
 
654
    if (pcdsc->choice == CddDescr_comment) {
 
655
      if (Nlm_StrCmp(pcdsc->data.ptrvalue,"linked to 3D-structure") != 0) {
 
656
        return (CharPtr) pcdsc->data.ptrvalue;
 
657
      }
 
658
    }
 
659
    pcdsc = pcdsc->next;
 
660
  }
 
661
  return NULL;
 
662
}
 
663
 
 
664
DatePtr LIBCALL CddGetCreateDate(CddPtr pcdd)
 
665
{
 
666
  CddDescrPtr pcdsc;
 
667
 
 
668
  pcdsc = pcdd->description;
 
669
  while (pcdsc) {
 
670
    if (pcdsc->choice == CddDescr_create_date) {
 
671
      return (DatePtr) pcdsc->data.ptrvalue;
 
672
    }
 
673
    pcdsc = pcdsc->next;
 
674
  }
 
675
  return NULL;
 
676
}
 
677
 
 
678
DatePtr LIBCALL CddGetUpdateDate(CddPtr pcdd)
 
679
{
 
680
  DatePtr pcd = NULL, pcdthis;
 
681
  CddDescrPtr pcdsc;
 
682
 
 
683
  pcdsc = pcdd->description;
 
684
  while (pcdsc) {
 
685
    if (pcdsc->choice == CddDescr_update_date) {
 
686
      pcdthis = (DatePtr) pcdsc->data.ptrvalue;
 
687
      if (!pcd) {
 
688
        pcd = pcdthis; 
 
689
      } else {
 
690
        if (pcdthis->data[1] > pcd->data[1]) {
 
691
          pcd = pcdthis;
 
692
        } else {
 
693
          if (pcdthis->data[2] > pcd->data[2]) {
 
694
            pcd = pcdthis;
 
695
          } else {
 
696
            if (pcdthis->data[3] > pcd->data[3]) pcd = pcdthis;
 
697
          }
 
698
        }
 
699
      }
 
700
    }
 
701
    pcdsc = pcdsc->next;
 
702
  }
 
703
  return pcd;
 
704
}
 
705
 
 
706
CharPtr LIBCALL CddGetSource(CddPtr pcdd)
 
707
{
 
708
  CddDescrPtr pcdsc;
 
709
 
 
710
  pcdsc = pcdd->description;
 
711
  while (pcdsc) {
 
712
    if (pcdsc->choice == CddDescr_source) {
 
713
      return (CharPtr) pcdsc->data.ptrvalue;
 
714
    }
 
715
    pcdsc = pcdsc->next;
 
716
  }
 
717
  return NULL;
 
718
}
 
719
 
 
720
CharPtr LIBCALL CddGetSourceId(CddPtr pcdd)
 
721
{
 
722
  CddDescrPtr pcdsc;
 
723
  GlobalIdPtr pGid;
 
724
  ValNodePtr  vnp;
 
725
 
 
726
  pcdsc = pcdd->description;
 
727
  while (pcdsc) {
 
728
    if (pcdsc->choice == CddDescr_source_id) {
 
729
      vnp = pcdsc->data.ptrvalue;
 
730
      if (vnp->choice == CddId_gid) {
 
731
        pGid = vnp->data.ptrvalue;
 
732
        return (CharPtr) pGid->accession;
 
733
      }
 
734
    }
 
735
    pcdsc = pcdsc->next;
 
736
  }
 
737
  return NULL;
 
738
}
 
739
 
 
740
Int4 LIBCALL CddGetAlignmentLength(CddPtr pcdd)
 
741
{
 
742
  SeqAnnotPtr sap;
 
743
  SeqAlignPtr salp;
 
744
  Int4        iLength = 0;
 
745
 
 
746
  sap = pcdd->seqannot;
 
747
  if (sap) {
 
748
    salp = (SeqAlignPtr) sap->data;
 
749
    while (salp) {
 
750
      iLength++;
 
751
      salp = salp->next;
 
752
    }
 
753
  }
 
754
  return iLength + 1;
 
755
}
 
756
 
 
757
/*---------------------------------------------------------------------------*/
 
758
/* return number of rows and average number of blocks per row                */
 
759
/*---------------------------------------------------------------------------*/
 
760
Int4Ptr LIBCALL GetAlignmentSize(SeqAlignPtr salp)
 
761
{
 
762
  Int4Ptr       pRet;
 
763
  Int4          iLength = 0;
 
764
  Int4          iBlocks = 0;
 
765
  DenseDiagPtr  ddp;
 
766
  SeqAlignPtr   salpThis;
 
767
 
 
768
  salpThis = salp;
 
769
  while (salpThis) {
 
770
    iLength++;
 
771
    ddp = salpThis->segs;
 
772
    while (ddp) {
 
773
      iBlocks++;
 
774
      ddp = ddp->next;
 
775
    }
 
776
    salpThis = salpThis->next;
 
777
  }
 
778
  pRet = MemNew(2*sizeof(Int4));
 
779
  pRet[0] = iLength + 1;
 
780
  pRet[1] = iBlocks / iLength;
 
781
  return pRet;
 
782
}
 
783
 
 
784
ValNodePtr LIBCALL CddGetAnnotNames(CddPtr pcdd)
 
785
{
 
786
  AlignAnnotPtr   aap;
 
787
  ValNodePtr      vnp, vnpHead = NULL;
 
788
 
 
789
  aap = pcdd->alignannot;
 
790
  while (aap) {
 
791
    vnp = ValNodeCopyStr(&(vnpHead),0,aap->description);
 
792
    aap = aap->next;  
 
793
  }
 
794
  return (vnpHead);
 
795
}
 
796
 
329
797
/*---------------------------------------------------------------------------*/
330
798
/*---------------------------------------------------------------------------*/
331
799
/* Query the status of a CD                                                  */
346
814
}
347
815
 
348
816
/*---------------------------------------------------------------------------*/
 
817
/* find a particular descriptive string in a CD                              */
 
818
/*---------------------------------------------------------------------------*/
 
819
Boolean LIBCALL CddHasDescription(CddPtr pcdd, CharPtr pc)
 
820
{
 
821
  CddDescrPtr    pCddesc;
 
822
  
 
823
  if (!pcdd) return FALSE;
 
824
  pCddesc = pcdd->description;
 
825
  while (pCddesc) {
 
826
    if (pCddesc->choice == CddDescr_comment) {
 
827
      if (Nlm_StrCmp((CharPtr)pCddesc->data.ptrvalue,pc) == 0) return TRUE;
 
828
    }
 
829
    pCddesc = pCddesc->next;
 
830
  }
 
831
  return FALSE;
 
832
}
 
833
 
 
834
/*---------------------------------------------------------------------------*/
 
835
/* Is the alignment decorated with feature annotation?                       */
 
836
/*---------------------------------------------------------------------------*/
 
837
Boolean LIBCALL CddHasAnnotation(CddPtr pcdd)
 
838
{
 
839
  AlignAnnotPtr aap;
 
840
  
 
841
  aap = pcdd->alignannot;
 
842
  if (aap) {
 
843
    if (aap->description) return TRUE;
 
844
  }
 
845
  return FALSE;
 
846
}
 
847
 
 
848
/*---------------------------------------------------------------------------*/
 
849
/*---------------------------------------------------------------------------*/
 
850
/* Is the CD master a 3D Structure                                           */
 
851
/*---------------------------------------------------------------------------*/
 
852
/*---------------------------------------------------------------------------*/
 
853
Boolean LIBCALL CddMasterIs3D(CddPtr pcdd)
 
854
{
 
855
  SeqAlignPtr     salp;
 
856
  DenseDiagPtr    ddp;
 
857
  SeqIdPtr        sip;
 
858
 
 
859
  if (!pcdd) return FALSE;
 
860
  if (!pcdd->seqannot) return FALSE;
 
861
  salp = pcdd->seqannot->data;
 
862
  ddp = salp->segs;
 
863
  sip = ddp->id;
 
864
  if (sip->choice == SEQID_PDB) return TRUE;
 
865
  if (sip->choice == SEQID_LOCAL) {
 
866
    sip = ddp->id->next;
 
867
    if (sip->choice == SEQID_PDB) return TRUE;
 
868
  }
 
869
  return FALSE;
 
870
}
 
871
 
 
872
/*---------------------------------------------------------------------------*/
 
873
/*---------------------------------------------------------------------------*/
 
874
/* How many structure-related alignments in a CDD?                           */
 
875
/*---------------------------------------------------------------------------*/
 
876
/*---------------------------------------------------------------------------*/
 
877
Int4 LIBCALL CddCount3DAlignments(CddPtr pcdd) 
 
878
{
 
879
  Int4              n3dali = 0;
 
880
  Boolean           bHasConsensus;
 
881
  Boolean           bHas3DMaster;
 
882
  SeqAlignPtr       salp;
 
883
  DenseDiagPtr      ddp;
 
884
  SeqIdPtr          sip;
 
885
  
 
886
  
 
887
  bHasConsensus = CddHasConsensus(pcdd);
 
888
  bHas3DMaster = CddMasterIs3D(pcdd);
 
889
 
 
890
  if (!bHas3DMaster) return(0);
 
891
  salp = pcdd->seqannot->data;
 
892
  while (salp) {
 
893
    ddp = salp->segs;
 
894
    sip = ddp->id->next;
 
895
    if (sip->choice == SEQID_PDB) n3dali++;
 
896
    salp = salp->next;
 
897
  }
 
898
  if (n3dali && bHasConsensus) n3dali--;
 
899
  return (n3dali);
 
900
}
 
901
 
 
902
/*---------------------------------------------------------------------------*/
349
903
/* Check if alignment master is a consensus Sequence                         */
350
904
/*---------------------------------------------------------------------------*/
351
905
Boolean LIBCALL SeqAlignHasConsensus(SeqAlignPtr salp)
367
921
/*---------------------------------------------------------------------------*/
368
922
/* Check if Cdd has a consensus Sequence                                     */
369
923
/*---------------------------------------------------------------------------*/
 
924
Boolean LIBCALL SipIsConsensus(SeqIdPtr sip)
 
925
{
 
926
  ObjectIdPtr   oidp;
 
927
  
 
928
  if (!sip) return (FALSE);
 
929
  if (sip->choice != SEQID_LOCAL) return (FALSE);
 
930
  oidp = sip->data.ptrvalue; if (!oidp) return(FALSE);
 
931
  if (StringCmp(oidp->str,"consensus")==0) return (TRUE);
 
932
  return(FALSE);
 
933
}
 
934
 
370
935
Boolean LIBCALL CddHasConsensus(CddPtr pcdd)
371
936
{
372
 
  SeqIntPtr   sintp;
373
 
  SeqIdPtr    sip;
374
 
  ObjectIdPtr oidp;
 
937
  SeqIntPtr     sintp;
 
938
  SeqIdPtr      sip;
 
939
  SeqAlignPtr   salp;
 
940
  DenseDiagPtr  ddp;
375
941
  
376
942
  if (!pcdd) return (FALSE);
377
943
  sintp = (SeqIntPtr) pcdd->profile_range;
378
 
  if (!sintp) return (FALSE);
 
944
  if (!sintp) {
 
945
    salp = pcdd->seqannot->data;
 
946
    ddp = salp->segs;
 
947
    sip = ddp->id;
 
948
    return (SipIsConsensus(sip));
 
949
  }
379
950
  sip = (SeqIdPtr) sintp->id;
380
 
  if (!sip) return (FALSE);
381
 
  if (sip->choice != SEQID_LOCAL) return (FALSE);
382
 
  oidp = sip->data.ptrvalue; if (!oidp) return(FALSE);
383
 
  if (StringCmp(oidp->str,"consensus")==0) return (TRUE);
384
 
  return(FALSE);
 
951
  return(SipIsConsensus(sip));
385
952
}
386
953
 
387
954
/*---------------------------------------------------------------------------*/
392
959
void LIBCALL CddRegularizeFileName(CharPtr cIn, CharPtr cAcc, CharPtr cFn,
393
960
                                   CharPtr cEx)
394
961
{
395
 
  CharPtr   cWhere;
396
 
  
397
962
  if (Nlm_StrStr(cIn, cEx) != 0) {
398
963
    Nlm_StrCpy(cFn,cIn);
399
964
    Nlm_StrNCpy(cAcc,cIn,Nlm_StrLen(cIn)-Nlm_StrLen(cEx));
413
978
 
414
979
/*---------------------------------------------------------------------------*/
415
980
/*---------------------------------------------------------------------------*/
 
981
/* repeat detection - find repeated subsequences in a CD consensus Sequence  */
 
982
/*---------------------------------------------------------------------------*/
 
983
/*---------------------------------------------------------------------------*/
 
984
Boolean LIBCALL CddCheckForRepeats(CddPtr pcdd, Int4 width, Int4 GapI, Int4 GapE,
 
985
                                   Nlm_FloatHi rthresh, Boolean bOutput)
 
986
{
 
987
  CddRepeatPtr        pcdr;
 
988
  BLAST_ScoreBlkPtr   sbp;
 
989
  BioseqPtr           bsp;
 
990
  Int4                len, winner = 1;
 
991
  Int4                **iscore, repct = 0;
 
992
  Int4                **fscore, bst, pnum, plen;
 
993
  Int2                iStatus;
 
994
  Int4                i, j, k, l, m, n, t1, t2, laststart;
 
995
  SeqPortPtr          spp;
 
996
  Uint1Ptr            buffer;
 
997
  Boolean             done = FALSE;
 
998
  Boolean             tbdone;
 
999
  Nlm_FloatHi         cfac, wplaus = 0.0;
 
1000
  SeqLocPtr           slp, slpHead;
 
1001
  SeqIntPtr           sintp;
 
1002
  
 
1003
  typedef struct repstep {
 
1004
    Int4            score;
 
1005
    Int4            start;
 
1006
    Int4            npred;
 
1007
    Int4            lpred;
 
1008
    Int4            lremn;
 
1009
    Nlm_FloatHi     plaus;
 
1010
  } RepStep;
 
1011
 
 
1012
  RepStep            rep[100];
 
1013
 
 
1014
  if (!CddHasConsensus) return FALSE;
 
1015
  bsp = pcdd->trunc_master;
 
1016
  len = bsp->length;
 
1017
  sbp = BLAST_ScoreBlkNew(Seq_code_ncbistdaa,1);
 
1018
  sbp->read_in_matrix = TRUE;
 
1019
  sbp->protein_alphabet = TRUE;
 
1020
  sbp->posMatrix = NULL;
 
1021
  sbp->number_of_contexts = 1;
 
1022
  iStatus = BlastScoreBlkMatFill(sbp,"BLOSUM62");
 
1023
  
 
1024
  spp = SeqPortNew(bsp, 0, LAST_RESIDUE, Seq_strand_unknown,Seq_code_ncbistdaa);
 
1025
  buffer = MemNew((len)*sizeof(Uint1));
 
1026
  for (i=0; i<len; i++) buffer[i] = SeqPortGetResidue(spp);
 
1027
  spp = SeqPortFree(spp);
 
1028
 
 
1029
  iscore = (Int4 **) MemNew(len * sizeof(Int4 *));
 
1030
  for (i=0;i<len;i++) iscore[i] = (Int4 *) MemNew(len * sizeof(Int4));
 
1031
  fscore = (Int4 **) MemNew(len * sizeof(Int4 *));
 
1032
  for (i=0;i<len;i++) fscore[i] = (Int4 *) MemNew(len * sizeof(Int4));
 
1033
 
 
1034
  for (i=0;i<len;i++) {
 
1035
    t1 = buffer[i];
 
1036
    iscore[i][i] = sbp->matrix[t1][t1];
 
1037
    for (j=i+1;j<len;j++) {
 
1038
      t2 = buffer[j];
 
1039
      iscore[i][j] = sbp->matrix[t1][t2];
 
1040
      fscore[j][i] = 0;
 
1041
    }  
 
1042
  }
 
1043
 
 
1044
  laststart = 0;
 
1045
  while (!done) {
 
1046
    for (i=0;i<len;i++) for (j=i;j<len;j++) fscore[i][j] = iscore[i][j];
 
1047
    for (i=len-2;i>=0;i--) {
 
1048
      for (j=len-2;j>=i;j--) {
 
1049
        bst = fscore[i+1][j+1];
 
1050
        for (k = i+2;k<i+2+width && k<len;k++) bst = MAX(bst,fscore[k][j+1]-GapI-GapE*(k-i-1));
 
1051
        for (l = j+2;l<j+2+width && l<len;l++) bst = MAX(bst,fscore[i+1][l]-GapI-GapE*(l-j-1));
 
1052
        bst = MAX(bst, 0);
 
1053
        fscore[i][j] += bst;
 
1054
      }
 
1055
    }
 
1056
 
 
1057
    bst = 0; l = 0;
 
1058
    for (j=laststart;j<len;j++) if (fscore[0][j] >= bst) {
 
1059
      bst = fscore[0][j]; l = j;
 
1060
    }
 
1061
    if (bst <= 0.0) {
 
1062
      done = TRUE;
 
1063
      break;
 
1064
    }
 
1065
    laststart = l;
 
1066
/*    printf("Next best alignment scores %4d at %4d (of %4d: %5.1f%%)\n",bst,l,len, 100.0*l/len); */
 
1067
    rep[repct].start = l;
 
1068
    rep[repct].score = bst;
 
1069
    rep[repct].npred = 1;
 
1070
    rep[repct].plaus = 0.0;
 
1071
    rep[repct].lpred = len;
 
1072
    rep[repct].lremn = len - rep[repct].start;
 
1073
    if (repct > 0) {
 
1074
      rep[repct].lpred = (Int4) (0.5 + ((Nlm_FloatHi)   l / (Nlm_FloatHi) repct));
 
1075
      rep[repct].npred = (Int4) (0.5 + ((Nlm_FloatHi) len / (Nlm_FloatHi) rep[repct].lpred)); 
 
1076
    }
 
1077
    repct++;
 
1078
    i = 0; j = l;
 
1079
    tbdone = FALSE;
 
1080
    if (l >= len-1) {
 
1081
      tbdone = TRUE;
 
1082
      done = TRUE;
 
1083
    }
 
1084
    while (!tbdone) {
 
1085
      m = i+1; n = j+1;
 
1086
      bst = fscore[m][n];
 
1087
      for (k=i+2;k<i+2+width && k<len;k++) {
 
1088
        if (fscore[k][j+1]-GapI-GapE*(k-i-1) > bst) {
 
1089
          bst = fscore[k][j+1]-GapI-GapE*(k-i-1);
 
1090
          m = k; n = j+1;
 
1091
        }
 
1092
      }
 
1093
      for (l=j+2;l<j+2+width && l<len;l++) {
 
1094
        if (fscore[i+1][l]-GapI-GapE*(l-j-1) > bst) {
 
1095
          bst = fscore[i+1][l]-GapI-GapE*(l-j-1);
 
1096
          m = i+1; n = l;
 
1097
        }
 
1098
      }
 
1099
      iscore[i][j] = -10000;
 
1100
      if (bst < 0) tbdone = TRUE;
 
1101
      i = m; j = n;
 
1102
      if (i >= len-1 || j >= len-1) {
 
1103
        tbdone = TRUE;
 
1104
        iscore[i][j] = -10000;
 
1105
      }
 
1106
    }
 
1107
  }
 
1108
  
 
1109
  for (i=0;i<len;i++) MemFree(fscore[i]);
 
1110
  MemFree(fscore);
 
1111
  for (i=0;i<len;i++) MemFree(iscore[i]);
 
1112
  MemFree(iscore);
 
1113
  rep[0].plaus = (Nlm_FloatHi) (rep[0].score - rep[1].score) / (Nlm_FloatHi) rep[0].score;
 
1114
  if (rep[0].plaus > rthresh) rep[0].plaus = 1.0;
 
1115
  winner = 1; wplaus = rep[0].plaus;
 
1116
 
 
1117
  for (i=1;i<repct;i++) {
 
1118
    pnum = i+1;
 
1119
    plen = (Int4) (((Nlm_FloatHi)len/(Nlm_FloatHi)pnum)+.5);
 
1120
    cfac = 0.0;
 
1121
    for (j=i;j>0;j--) {
 
1122
      cfac += (Nlm_FloatHi)(abs(rep[j].start-(j*plen)))/(Nlm_FloatHi)plen;
 
1123
    }
 
1124
    cfac /= (Nlm_FloatHi) i;
 
1125
    rep[i].plaus = 1.0 - cfac;
 
1126
 
 
1127
    if (rep[i].plaus > wplaus) {
 
1128
      wplaus = rep[i].plaus;
 
1129
      winner = pnum;
 
1130
    }
 
1131
  }
 
1132
  for (i=0;i<winner;i++) {
 
1133
    if (rep[i].lpred < 3 || (i>0 && rep[i].score < 17)) {
 
1134
      winner = 1; break;
 
1135
    }
 
1136
  }
 
1137
 
 
1138
  if (bOutput) {
 
1139
    for (i=0;i<repct;i++) {
 
1140
      printf("%s %3d %4d %8.5f %4d %8.5f %2d %3d %8.5f\n",pcdd->name,
 
1141
             i,rep[i].score,(Nlm_FloatHi)rep[i].score/(Nlm_FloatHi)rep[0].score,
 
1142
             rep[i].start, (Nlm_FloatHi)rep[i].start/(Nlm_FloatHi)len,
 
1143
             rep[i].npred,rep[i].lpred, rep[i].plaus);
 
1144
    }
 
1145
    printf("%s: predicted %d repeats\n",pcdd->name,winner);
 
1146
  }
 
1147
 
 
1148
  CddKillDescr(pcdd,NULL,CddDescr_repeats,0);
 
1149
  if (winner > 1) {
 
1150
    pcdr = CddRepeatNew(); pcdr->avglen = 0;  
 
1151
    slp = NULL; slpHead = NULL;
 
1152
    for (i=0;i<winner;i++) {
 
1153
      k = rep[i].start;
 
1154
      if (i < (winner-1)) l = rep[i+1].start - 1; else l = len - 1;
 
1155
      pcdr->avglen += l-k+1;
 
1156
      slp = (SeqLocPtr) ValNodeNew(NULL);
 
1157
      slp->choice = SEQLOC_INT;
 
1158
      sintp = SeqIntNew();
 
1159
      sintp->from = k;
 
1160
      sintp->to = l;
 
1161
      sintp->id = SeqIdDup(pcdd->profile_range->id);
 
1162
      slp->data.ptrvalue = sintp;
 
1163
      if (!slpHead) slpHead = slp; else ValNodeLink(&(slpHead),slp);
 
1164
    }
 
1165
    pcdr->avglen /= winner;
 
1166
    pcdr->count = winner;
 
1167
    pcdr->location = (SeqLocPtr) ValNodeNew(NULL);
 
1168
    pcdr->location->choice = SEQLOC_PACKED_INT;
 
1169
    pcdr->location->data.ptrvalue = slpHead;
 
1170
    CddAssignDescr(pcdd,pcdr,CddDescr_repeats,0);
 
1171
  }
 
1172
  
 
1173
  sbp = BLAST_ScoreBlkDestruct(sbp);
 
1174
  MemFree(buffer);
 
1175
  if (winner > 1) return TRUE;
 
1176
  return (FALSE);
 
1177
}
 
1178
 
 
1179
 
 
1180
/*---------------------------------------------------------------------------*/
 
1181
/*---------------------------------------------------------------------------*/
416
1182
/* retrieve the accession of a CD as a character string                      */
417
1183
/*---------------------------------------------------------------------------*/
418
1184
/*---------------------------------------------------------------------------*/
439
1205
/* report Errors in processing and exit immediately                          */
440
1206
/*---------------------------------------------------------------------------*/
441
1207
/*---------------------------------------------------------------------------*/
442
 
void LIBCALL CddHtmlError(CharPtr cErrTxt) 
 
1208
void LIBCALL CddSimpleHtmlError(CharPtr cErrTxt) 
443
1209
{
444
1210
  printf("Content-type: text/html\n\n");
445
1211
  printf("<h2>Error:</h2>\n");
493
1259
  Uint1       seqtype;
494
1260
  ValNodePtr  tmp;
495
1261
  ObjectIdPtr oid;
496
 
  Int4        len, i;
497
 
  Int2        residue;
 
1262
  Int4        len;
498
1263
  ValNode     fake;
499
1264
  SeqLocPtr   the_segs, head, curr;
500
1265
  Boolean     handled = FALSE, split;
537
1302
 
538
1303
  if ((newbsp->repr == Seq_repr_raw) ||
539
1304
    (newbsp->repr == Seq_repr_const)) {
540
 
      AsnTypePtr  atp;
541
 
      AsnIoPtr  aip;
542
1305
    if (ISA_aa(newbsp->mol)) {
543
1306
      seqtype = Seq_code_ncbieaa;
544
1307
    } else {
656
1419
 
657
1420
/*---------------------------------------------------------------------------*/
658
1421
/*---------------------------------------------------------------------------*/
 
1422
/* Create a Dense-Diag seqalign from a Discontinuous SeqAlign                */
 
1423
/*---------------------------------------------------------------------------*/
 
1424
/*---------------------------------------------------------------------------*/
 
1425
SeqAlignPtr LIBCALL CddMSLMixedToMSLDenDiag(SeqAlignPtr salp)
 
1426
{
 
1427
  SeqAlignPtr         salpSub, salpnew, salphead, salptail;
 
1428
  DenseDiagPtr        ddp, ddphead, ddptail;
 
1429
  DenseSegPtr         dsp;
 
1430
  Int4                i, s1, s2;
 
1431
  ScorePtr            scp, scpHead, scpTail, scpThis;
 
1432
 
 
1433
  if (!salp) return(NULL);
 
1434
  salphead = NULL;
 
1435
  salptail = NULL;
 
1436
  while (salp) {
 
1437
    if (salp->segtype) {
 
1438
      if (salp->segtype != SAS_DISC && salp->segtype != SAS_DENSEG) {
 
1439
        CddSevError("CddMSLDiscToMSLDenDiag: Wrong segtype specified!");
 
1440
      }
 
1441
    }
 
1442
    if (salp->dim) {
 
1443
      if (salp->dim != 2) {
 
1444
        CddSevError("CddMSLDiscToMSLDenDiag: Expect alignments of dimension 2!");
 
1445
      }
 
1446
    }
 
1447
    ddphead = NULL;
 
1448
    ddptail = NULL;
 
1449
 
 
1450
    if (salp->segtype == SAS_DISC) {
 
1451
      salpSub = salp->segs;
 
1452
      while (salpSub) {
 
1453
        if (salpSub->segtype != SAS_DENSEG)
 
1454
          CddSevError("CddMSLDiscToMSLDenDiag: Wrong segtype in Sub-Alignment!");
 
1455
        if (salpSub->dim && salpSub->dim !=2)
 
1456
          CddSevError("CddMSLDiscToMSLDenDiag: Wrong dimension in Sub-Alignment!");
 
1457
        dsp = salpSub->segs;
 
1458
        for (i=0;i<dsp->numseg;i++) {
 
1459
          s1 = dsp->starts[i*2];
 
1460
          s2 = dsp->starts[i*2+1]; 
 
1461
          if (s1 >=0 && s2>= 0) {
 
1462
            ddp = DenseDiagNew();
 
1463
            ddp->starts = MemNew(2*sizeof(Int4));
 
1464
            ddp->starts[0] = s1;
 
1465
            ddp->starts[1] = s2;
 
1466
            ddp->len=dsp->lens[i];
 
1467
            ddp->id = SeqIdDupList(dsp->ids);
 
1468
            ddp->dim = 2;
 
1469
            if (!ddphead) {
 
1470
              ddphead = ddp;
 
1471
              ddptail = ddp;
 
1472
              ddptail->next = NULL;
 
1473
            } else {
 
1474
              ddptail->next = ddp;
 
1475
              ddptail = ddp;
 
1476
              ddptail->next = NULL;
 
1477
            }
 
1478
          }
 
1479
        }
 
1480
        salpSub = salpSub->next;  
 
1481
      }
 
1482
    } else { /* segtype is plain SAS_DENSEG */
 
1483
      dsp = salp->segs;
 
1484
      ddphead = NULL;
 
1485
      ddptail = NULL;
 
1486
      for (i=0;i<dsp->numseg;i++) {
 
1487
        s1 = dsp->starts[i*2];
 
1488
        s2 = dsp->starts[i*2+1]; 
 
1489
        if (s1 >=0 && s2>= 0) {
 
1490
          ddp = DenseDiagNew();
 
1491
          ddp->starts = MemNew(2*sizeof(Int4));
 
1492
          ddp->starts[0] = s1;
 
1493
          ddp->starts[1] = s2;
 
1494
          ddp->len=dsp->lens[i];
 
1495
          ddp->id = SeqIdDupList(dsp->ids);
 
1496
          ddp->dim = 2;
 
1497
          if (!ddphead) {
 
1498
            ddphead = ddp;
 
1499
            ddptail = ddp;
 
1500
            ddptail->next = NULL;
 
1501
          } else {
 
1502
            ddptail->next = ddp;
 
1503
            ddptail = ddp;
 
1504
            ddptail->next = NULL;
 
1505
          }
 
1506
        }
 
1507
      }
 
1508
    }
 
1509
    salpnew = SeqAlignNew();
 
1510
    salpnew->type = SAT_PARTIAL;
 
1511
    salpnew->segtype = SAS_DENDIAG;
 
1512
    salpnew->dim = 2;
 
1513
    salpnew->segs = ddphead;
 
1514
    scp = salp->score; scpHead = NULL;
 
1515
    while (scp) {
 
1516
      scpThis = ScoreNew();
 
1517
      scpThis->id = ObjectIdDup(scp->id);
 
1518
      scpThis->choice = scp->choice;
 
1519
      scpThis->value = scp->value;
 
1520
      if (!scpHead) {
 
1521
        scpHead = scpThis;
 
1522
      } else {
 
1523
        scpTail->next = scpThis;
 
1524
      }
 
1525
      scpTail = scpThis;
 
1526
      scp = scp->next;
 
1527
    }
 
1528
    salpnew->score = scpHead;
 
1529
    if (!salphead) {
 
1530
      salphead = salpnew;
 
1531
      salptail = salpnew;
 
1532
      salptail->next = NULL;
 
1533
    } else {
 
1534
      salptail->next = salpnew;
 
1535
      salptail = salpnew;
 
1536
      salptail->next = NULL;
 
1537
    }
 
1538
    salp = salp->next;
 
1539
  }
 
1540
  return(salphead);
 
1541
}
 
1542
 
 
1543
/*---------------------------------------------------------------------------*/
 
1544
/*---------------------------------------------------------------------------*/
659
1545
/* Create a Dense-Diag seqalign from a Dense-Seg SeqAlign                    */
660
1546
/*---------------------------------------------------------------------------*/
661
1547
/*---------------------------------------------------------------------------*/
665
1551
  DenseDiagPtr        ddp, ddphead, ddptail;
666
1552
  DenseSegPtr         dsp;
667
1553
  Int4                i, s1, s2;
 
1554
  ScorePtr            scp, scpHead, scpTail, scpThis;
668
1555
 
669
1556
 
670
1557
  if (!salp) return(NULL);
693
1580
        ddp->starts[0] = s1;
694
1581
        ddp->starts[1] = s2;
695
1582
        ddp->len=dsp->lens[i];
696
 
        ddp->id = dsp->ids;
 
1583
        ddp->id = SeqIdDupList(dsp->ids);
697
1584
        ddp->dim = 2;
698
1585
        if (!ddphead) {
699
1586
          ddphead = ddp;
711
1598
    salpnew->segtype = SAS_DENDIAG;
712
1599
    salpnew->dim = 2;
713
1600
    salpnew->segs = ddphead;
 
1601
    scp = salp->score; scpHead = NULL;
 
1602
    while (scp) {
 
1603
      scpThis = ScoreNew();
 
1604
      scpThis->id = ObjectIdDup(scp->id);
 
1605
      scpThis->choice = scp->choice;
 
1606
      scpThis->value = scp->value;
 
1607
      if (!scpHead) {
 
1608
        scpHead = scpThis;
 
1609
      } else {
 
1610
        scpTail->next = scpThis;
 
1611
      }
 
1612
      scpTail = scpThis;
 
1613
      scp = scp->next;
 
1614
    }
 
1615
    salpnew->score = scpHead;
714
1616
    if (!salphead) {
715
1617
      salphead = salpnew;
716
1618
      salptail = salpnew;
726
1628
}
727
1629
 
728
1630
/*---------------------------------------------------------------------------*/
 
1631
/*---------------------------------------------------------------------------*/
729
1632
/* make a DenseSeg Alignment-Set turning each diag into a separate alignment */
730
1633
/*---------------------------------------------------------------------------*/
 
1634
/*---------------------------------------------------------------------------*/
731
1635
SeqAlignPtr LIBCALL CddMSLDenDiagToMSLDenSeg(SeqAlignPtr salp)
732
1636
{
733
1637
  SeqAlignPtr         salpnew, salphead, salptail;
763
1667
      salpnew->segs = dsp;
764
1668
      dsp->dim = 2;
765
1669
      dsp->numseg = 1;
766
 
      dsp->ids = ddp->id;
 
1670
      dsp->ids = SeqIdDupList(ddp->id);
767
1671
      dsp->starts = ddp->starts;
768
1672
      dsp->lens = &(ddp->len);
769
1673
      dsp->strands = ddp->strands;
778
1682
}
779
1683
 
780
1684
/*---------------------------------------------------------------------------*/
 
1685
/*---------------------------------------------------------------------------*/
781
1686
/* convert a dendiag pairwise alignment set into a multiple dendiag alignment*/
782
1687
/* if possible - i.e. check that the number of segments and their extents on */
783
1688
/* the common master is the same for all pairwise alignments                 */
784
1689
/*---------------------------------------------------------------------------*/
 
1690
/*---------------------------------------------------------------------------*/
785
1691
SeqAlignPtr LIBCALL CddMSLDenDiagToMULDenDiag(SeqAlignPtr salp)
786
1692
{
787
1693
  SeqAlignPtr     salpHead;
910
1816
/* Calculate per column information content for a SeqAlign                   */
911
1817
/*---------------------------------------------------------------------------*/
912
1818
/*---------------------------------------------------------------------------*/
 
1819
Nlm_FloatHi LIBCALL SeqAlignConservation(SeqAlignPtr salp, Nlm_FloatHi fract,
 
1820
                                         BioseqPtr bsp_master, Boolean bHasConsensus, Int4 offset)
 
1821
{
 
1822
  Int4               i, c, a;
 
1823
  Int4               qlength, nover = 0;
 
1824
  Int4               alphabetSize = 26;
 
1825
  Int4               nColumns;
 
1826
  Int4               *ntypes;
 
1827
  Nlm_FloatHi        **typefreq, sumfreq;
 
1828
  DenseDiagPtr       ddp;
 
1829
  BioseqPtr          bsp;
 
1830
  SeqIdPtr           sip;
 
1831
  SeqPortPtr         spp;
 
1832
  Uint1Ptr           buffer;
 
1833
  
 
1834
  if (!salp) return(0.0);
 
1835
  if (!bsp_master) return(0.0);
 
1836
  qlength = bsp_master->length;
 
1837
  ntypes = (Int4 *) MemNew(qlength*sizeof(Int4));
 
1838
  for (i=0;i<qlength;i++) ntypes[i] = 0;
 
1839
  typefreq = (Nlm_FloatHi**) MemNew(qlength * sizeof(Nlm_FloatHi *));
 
1840
  for (i=0;i<qlength;i++) {
 
1841
    typefreq[i] = (Nlm_FloatHi *)MemNew(alphabetSize * sizeof(Nlm_FloatHi));
 
1842
    for (a=0;a<alphabetSize;a++) typefreq[i][a] = 0.0;
 
1843
  }
 
1844
  if (!bHasConsensus) {   /* count residues in the master/representative too */
 
1845
    ddp = salp->segs;
 
1846
    sip = ddp->id;
 
1847
    bsp = bsp_master;
 
1848
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
1849
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
1850
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
1851
    spp = SeqPortFree(spp);
 
1852
    while (ddp) {
 
1853
      for (c=ddp->starts[0]-offset;c<ddp->starts[0]-offset+ddp->len;c++) {
 
1854
        if (buffer[c] >= 0 && buffer[c] < (Uint1) alphabetSize)
 
1855
          typefreq[c][buffer[c]] += 1.0;
 
1856
      }
 
1857
      ddp = ddp->next;
 
1858
    }
 
1859
    MemFree(buffer);
 
1860
  }  
 
1861
  while (salp) {
 
1862
    ddp = salp->segs;
 
1863
    sip = ddp->id->next;
 
1864
    bsp = BioseqLockById(sip);
 
1865
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
1866
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
1867
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
1868
    spp = SeqPortFree(spp);
 
1869
    BioseqUnlock(bsp);
 
1870
    while (ddp) {
 
1871
      for (c=ddp->starts[1];c<ddp->starts[1]+ddp->len;c++) {
 
1872
        i = ddp->starts[0]-offset+c-ddp->starts[1];
 
1873
        if (buffer[c] >= 0 && buffer[c] < (Uint1) alphabetSize)
 
1874
          typefreq[i][buffer[c]] += 1.0;
 
1875
      }
 
1876
      ddp = ddp->next;
 
1877
    }
 
1878
    MemFree(buffer);
 
1879
    salp = salp->next;
 
1880
  }
 
1881
  for (i=0;i<qlength;i++) {
 
1882
    sumfreq = 0.0;
 
1883
    for (a=0;a<alphabetSize;a++) {
 
1884
      if (typefreq[i][a] > 0.0) {
 
1885
        ntypes[i]++;
 
1886
        sumfreq += typefreq[i][a];
 
1887
      }
 
1888
    }
 
1889
    if (sumfreq > 0.0) for (a=0;a<alphabetSize;a++) {
 
1890
      typefreq[i][a] /= sumfreq;
 
1891
    }
 
1892
    for (a=0;a<alphabetSize;a++) {
 
1893
      if (typefreq[i][a] >= fract) {
 
1894
        nover++; break;
 
1895
      } 
 
1896
    }
 
1897
  }
 
1898
  for (i=0;i<qlength;i++) MemFree(typefreq[i]);
 
1899
  MemFree(typefreq);
 
1900
  MemFree(ntypes);
 
1901
  return((Nlm_FloatHi)nover/(Nlm_FloatHi)qlength);
 
1902
}
 
1903
 
 
1904
/*---------------------------------------------------------------------------*/
 
1905
/*---------------------------------------------------------------------------*/
 
1906
/* Calculate per column information content for a SeqAlign                   */
 
1907
/*---------------------------------------------------------------------------*/
 
1908
/*---------------------------------------------------------------------------*/
913
1909
Nlm_FloatHiPtr LIBCALL SeqAlignInform(SeqAlignPtr salp, BioseqPtr bsp_master,
914
 
                                      Boolean bHasConsensus)
 
1910
                                      Boolean bHasConsensus, Int4 offset)
915
1911
{
916
1912
  BLAST_ScoreBlkPtr  sbp;
917
1913
  BLAST_ResFreqPtr   stdrfp;
921
1917
  Int4               i, c, a;
922
1918
  Int4               qlength;
923
1919
  Int4               alphabetSize = 26;
924
 
  MatrixPtr          posfreq;
925
 
  ValNodePtr         vnp;
926
 
  Int4               nColumns;
927
1920
  Nlm_FloatHi        posEpsilon = 0.0001;
928
1921
  Nlm_FloatHi        QoverPest;
929
1922
  Int4               *ntypes;
933
1926
  SeqIdPtr           sip;
934
1927
  SeqPortPtr         spp;
935
1928
  Uint1Ptr           buffer;
 
1929
  Char               matrix[32];
936
1930
  
937
1931
  if (!salp) return(NULL);
938
1932
  if (!bsp_master) return(NULL);
943
1937
  sbp->protein_alphabet = TRUE;
944
1938
  sbp->posMatrix = NULL;
945
1939
  sbp->number_of_contexts = 1;
946
 
  iStatus = BlastScoreBlkMatFill(sbp,"BLOSUM62");
 
1940
  StrCpy(matrix,"BLOSUM62");
 
1941
  iStatus = BlastScoreBlkMatFill(sbp,matrix);
947
1942
  stdrfp = BlastResFreqNew(sbp);
948
1943
  BlastResFreqStdComp(sbp,stdrfp); 
949
1944
  standardProb = MemNew(alphabetSize * sizeof(Nlm_FloatHi));
950
1945
  for(a = 0; a < alphabetSize; a++) standardProb[a] = stdrfp->prob[a];
951
1946
  stdrfp = BlastResFreqDestruct(stdrfp);
952
 
  ntypes = (Int4 *) MemNew(qlength*sizeof(Int4));
953
 
  for (i=0;i<qlength;i++) ntypes[i] = 0;
954
 
  typefreq = (Nlm_FloatHi**) MemNew(qlength * sizeof(Nlm_FloatHi *));
955
 
  for (i=0;i<qlength;i++) {
956
 
    typefreq[i] = (Nlm_FloatHi *)MemNew(alphabetSize * sizeof(Nlm_FloatHi));
957
 
    for (a=0;a<alphabetSize;a++) typefreq[i][a] = 0.0;
958
 
  }
959
 
  if (!bHasConsensus) {   /* count residues in the master/representative too */
960
 
    ddp = salp->segs;
961
 
    sip = ddp->id;
962
 
    bsp = bsp_master;
963
 
    buffer = MemNew((bsp->length)*sizeof(Uint1));
964
 
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
965
 
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
966
 
    spp = SeqPortFree(spp);
967
 
    while (ddp) {
968
 
      for (c=ddp->starts[0];c<ddp->starts[0]+ddp->len;c++) {
969
 
        if (buffer[c] >= 0 && buffer[c] < alphabetSize)
970
 
          typefreq[c][buffer[c]] += 1.0;
971
 
      }
972
 
      ddp = ddp->next;
973
 
    }
974
 
    MemFree(buffer);
975
 
  }  
976
 
  while (salp) {
977
 
    ddp = salp->segs;
978
 
    sip = ddp->id->next;
979
 
    bsp = CddRetrieveBioseqById(sip,NULL);
980
 
    buffer = MemNew((bsp->length)*sizeof(Uint1));
981
 
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
982
 
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
983
 
    spp = SeqPortFree(spp);
984
 
    while (ddp) {
985
 
      for (c=ddp->starts[1];c<ddp->starts[1]+ddp->len;c++) {
986
 
        i = ddp->starts[0]+c-ddp->starts[1];
987
 
        if (buffer[c] >= 0 && buffer[c] < alphabetSize)
988
 
          typefreq[i][buffer[c]] += 1.0;
989
 
      }
990
 
      ddp = ddp->next;
991
 
    }
992
 
    MemFree(buffer);
993
 
    salp = salp->next;
994
 
  }
995
 
  for (i=0;i<qlength;i++) {
996
 
    sumfreq = 0.0;
997
 
    for (a=0;a<alphabetSize;a++) {
998
 
      if (typefreq[i][a] > 0.0) {
999
 
        ntypes[i]++;
1000
 
        sumfreq += typefreq[i][a];
1001
 
      }
1002
 
    }
1003
 
    if (sumfreq > 0.0) for (a=0;a<alphabetSize;a++) {
1004
 
      typefreq[i][a] /= sumfreq; 
1005
 
    }
1006
 
  }
1007
 
 
1008
 
  for (c=0;c<qlength;c++) {
1009
 
    Informativeness[c] = 0.0;
1010
 
    for (a=0;a<alphabetSize;a++) {
1011
 
      if (standardProb[a] > posEpsilon) {
1012
 
        QoverPest = typefreq[c][a] / standardProb[a];
1013
 
        if (QoverPest > posEpsilon) {
1014
 
          Informativeness[c] += typefreq[c][a] * log(QoverPest) / NCBIMATH_LN2;        
1015
 
        }
1016
 
      }
1017
 
    }
1018
 
  }
1019
 
  for (i=0;i<qlength;i++) MemFree(typefreq[i]);
1020
 
  MemFree(typefreq);
1021
 
  MemFree(ntypes);
1022
 
  MemFree(standardProb);
1023
1947
  sbp = (BLAST_ScoreBlkPtr) BLAST_ScoreBlkDestruct(sbp);
 
1948
  ntypes = (Int4 *) MemNew(qlength*sizeof(Int4));
 
1949
  for (i=0;i<qlength;i++) ntypes[i] = 0;
 
1950
  typefreq = (Nlm_FloatHi**) MemNew(qlength * sizeof(Nlm_FloatHi *));
 
1951
  for (i=0;i<qlength;i++) {
 
1952
    typefreq[i] = (Nlm_FloatHi *)MemNew(alphabetSize * sizeof(Nlm_FloatHi));
 
1953
    for (a=0;a<alphabetSize;a++) typefreq[i][a] = 0.0;
 
1954
  }
 
1955
  if (!bHasConsensus) {   /* count residues in the master/representative too */
 
1956
    ddp = salp->segs;
 
1957
    sip = ddp->id;
 
1958
    bsp = bsp_master;
 
1959
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
1960
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
1961
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
1962
    spp = SeqPortFree(spp);
 
1963
    while (ddp) {
 
1964
      for (c=ddp->starts[0]-offset;c<ddp->starts[0]-offset+ddp->len;c++) {
 
1965
        if (buffer[c] < (Uint1) alphabetSize)
 
1966
          typefreq[c][buffer[c]] += 1.0;
 
1967
      }
 
1968
      ddp = ddp->next;
 
1969
    }
 
1970
    MemFree(buffer);
 
1971
  }  
 
1972
  while (salp) {
 
1973
    ddp = salp->segs;
 
1974
    sip = ddp->id->next;
 
1975
    bsp = BioseqLockById(sip);
 
1976
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
1977
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
1978
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
1979
    spp = SeqPortFree(spp);
 
1980
    BioseqUnlock(bsp);
 
1981
    while (ddp) {
 
1982
      for (c=ddp->starts[1];c<ddp->starts[1]+ddp->len;c++) {
 
1983
        i = ddp->starts[0]-offset+c-ddp->starts[1];
 
1984
        if (buffer[c] < (Uint1) alphabetSize)
 
1985
          typefreq[i][buffer[c]] += 1.0;
 
1986
      }
 
1987
      ddp = ddp->next;
 
1988
    }
 
1989
    MemFree(buffer);
 
1990
    salp = salp->next;
 
1991
  }
 
1992
  for (i=0;i<qlength;i++) {
 
1993
    sumfreq = 0.0;
 
1994
    for (a=0;a<alphabetSize;a++) {
 
1995
      if (typefreq[i][a] > 0.0) {
 
1996
        ntypes[i]++;
 
1997
        sumfreq += typefreq[i][a];
 
1998
      }
 
1999
    }
 
2000
    if (sumfreq > 0.0) for (a=0;a<alphabetSize;a++) {
 
2001
      typefreq[i][a] /= sumfreq; 
 
2002
    }
 
2003
  }
 
2004
 
 
2005
  for (c=0;c<qlength;c++) {
 
2006
    Informativeness[c] = 0.0;
 
2007
    for (a=0;a<alphabetSize;a++) {
 
2008
      if (standardProb[a] > posEpsilon) {
 
2009
        QoverPest = typefreq[c][a] / standardProb[a];
 
2010
        if (QoverPest > posEpsilon) {
 
2011
          Informativeness[c] += typefreq[c][a] * log(QoverPest) / NCBIMATH_LN2;        
 
2012
        }
 
2013
      }
 
2014
    }
 
2015
  }
 
2016
  for (i=0;i<qlength;i++) MemFree(typefreq[i]);
 
2017
  MemFree(typefreq);
 
2018
  MemFree(ntypes);
 
2019
  MemFree(standardProb);
1024
2020
  return(Informativeness);
1025
2021
}
1026
2022
 
1029
2025
/* Calculate per column information content for a CDD from the seqalign      */
1030
2026
/*---------------------------------------------------------------------------*/
1031
2027
/*---------------------------------------------------------------------------*/
 
2028
Int4 ** LIBCALL CddCountResTypes(CddPtr pcdd, Int4 *ncols)
 
2029
{
 
2030
  BioseqSetPtr       bssp;
 
2031
  Boolean            bHasConsensus = CddHasConsensus(pcdd);
 
2032
  Int4               offset, qlength, i, j, c;
 
2033
  Int4             **ResFreqTable;
 
2034
  SeqAlignPtr        salp;
 
2035
  DenseDiagPtr       ddp;
 
2036
  SeqIdPtr           sip;
 
2037
  BioseqPtr          bsp;
 
2038
  Uint1             *buffer;
 
2039
  SeqPortPtr         spp;
 
2040
 
 
2041
  if (!pcdd) return(NULL);
 
2042
  if (!pcdd->trunc_master) return(NULL);
 
2043
  bssp = pcdd->sequences->data.ptrvalue;
 
2044
  offset = pcdd->profile_range->from;
 
2045
  qlength = pcdd->trunc_master->length;
 
2046
  *ncols = qlength;  
 
2047
  ResFreqTable = MemNew(qlength*sizeof(Int4 *));
 
2048
  for (i=0;i<qlength;i++) {
 
2049
    ResFreqTable[i] = MemNew(26*sizeof(Int4));
 
2050
    for (j=0;j<26;j++) ResFreqTable[i][j] = 0;
 
2051
  }
 
2052
  salp = pcdd->seqannot->data;
 
2053
  if (!bHasConsensus) {   /* count residues in the master/representative too */
 
2054
    ddp = salp->segs;
 
2055
    sip = ddp->id;
 
2056
    bsp = CddRetrieveBioseqById(sip,bssp->seq_set);
 
2057
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
2058
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
2059
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
2060
    spp = SeqPortFree(spp);
 
2061
    while (ddp) {
 
2062
      for (c=ddp->starts[0];c<ddp->starts[0]+ddp->len;c++) {
 
2063
        ResFreqTable[c-offset][buffer[c]]++;
 
2064
      }
 
2065
      ddp = ddp->next;
 
2066
    }
 
2067
    MemFree(buffer);
 
2068
  }  
 
2069
  while (salp) {
 
2070
    ddp = salp->segs;
 
2071
    sip = ddp->id->next;
 
2072
    bsp = CddRetrieveBioseqById(sip,bssp->seq_set);
 
2073
    buffer = MemNew((bsp->length)*sizeof(Uint1));
 
2074
    spp = SeqPortNew(bsp, 0, bsp->length-1, Seq_strand_unknown, Seq_code_ncbistdaa);
 
2075
    for (i=0; i<bsp->length;i++) buffer[i] = SeqPortGetResidue(spp);
 
2076
    spp = SeqPortFree(spp);
 
2077
    while (ddp) {
 
2078
      for (c=ddp->starts[1];c<ddp->starts[1]+ddp->len;c++) {
 
2079
        i = ddp->starts[0]-offset+c-ddp->starts[1];
 
2080
        ResFreqTable[i][buffer[c]]++;
 
2081
      }
 
2082
      ddp = ddp->next;
 
2083
    }
 
2084
    MemFree(buffer);
 
2085
    salp = salp->next;
 
2086
  }
 
2087
  return(ResFreqTable);
 
2088
}
 
2089
 
 
2090
/*---------------------------------------------------------------------------*/
 
2091
/*---------------------------------------------------------------------------*/
 
2092
/* Calculate per column information content for a CDD from the seqalign      */
 
2093
/*---------------------------------------------------------------------------*/
 
2094
/*---------------------------------------------------------------------------*/
1032
2095
Nlm_FloatHiPtr LIBCALL CddAlignInform(CddPtr pcdd, Nlm_FloatHi * Niobs)
1033
2096
{
1034
2097
  Int4               offset;
1041
2104
  Int4               i, c, a;
1042
2105
  Int4               qlength;
1043
2106
  Int4               alphabetSize = 26;
1044
 
  MatrixPtr          posfreq;
1045
 
  ValNodePtr         vnp;
1046
 
  Int4               nColumns;
1047
2107
  Nlm_FloatHi        posEpsilon = 0.0001;
1048
2108
  Nlm_FloatHi        QoverPest;
1049
2109
  Int4               *ntypes;
1164
2224
  Nlm_FloatHiPtr     Informativeness;
1165
2225
  Nlm_FloatHi        thisposFreq;
1166
2226
  Int2               iStatus;
1167
 
  Int4               i, c, a;
 
2227
  Int4               c, a;
1168
2228
  Int4               qlength;
1169
2229
  Int4               alphabetSize = 26;
1170
2230
  MatrixPtr          posfreq;
1171
2231
  ValNodePtr         vnp;
1172
 
  Int4               nColumns;
1173
2232
  Nlm_FloatHi        posEpsilon = 0.0001;
1174
2233
  Nlm_FloatHi        QoverPest;
1175
2234
 
1223
2282
  Nlm_FloatHiPtr     Informativeness;
1224
2283
  Nlm_FloatHi        thisposFreq;
1225
2284
  Int2               iStatus;
1226
 
  Int4               i, c, a;
 
2285
  Int4               c, a;
1227
2286
  Int4               qlength;
1228
2287
  Int4               alphabetSize = 26;
1229
 
  Int4               nColumns;
1230
2288
  Nlm_FloatHi        posEpsilon = 0.0001;
1231
2289
  Nlm_FloatHi        QoverPest;
1232
2290
 
1778
2836
 
1779
2837
  while (sepThis) {
1780
2838
    if (IS_Bioseq(sepThis)) {
1781
 
      bsp = sep->data.ptrvalue;
 
2839
      bsp = sepThis->data.ptrvalue;
1782
2840
      if (CddSameSip(bsp->id, sip)) return(bsp);
1783
2841
    } else if (IS_Bioseq_set(sepThis)) {
1784
 
      bssp = sep->data.ptrvalue;
 
2842
      bssp = sepThis->data.ptrvalue;
1785
2843
      bsp = CddFindBioseqInSeqEntry(bssp->seq_set,sip);
1786
2844
      if (CddSameSip(bsp->id, sip)) return(bsp);
1787
2845
    }
1968
3026
  pCDea->length = 0;
1969
3027
  pCDea->ids = NULL;
1970
3028
  pCDea->adata = NULL;
 
3029
  pCDea->starts = NULL;
1971
3030
  pCDea->bIdAlloc = FALSE;
1972
3031
  return(pCDea);
1973
3032
}
1974
3033
 
1975
3034
CddExpAlignPtr CddExpAlignFree(CddExpAlignPtr pCDea)
1976
3035
{
1977
 
  SeqIdPtr sip, sipNext;
1978
 
 
1979
3036
  if (!pCDea) return NULL;
1980
3037
  if (NULL != pCDea->adata) free(pCDea->adata);
1981
3038
  pCDea->ids = SeqIdSetFree(pCDea->ids);
 
3039
  if (NULL != pCDea->starts) free(pCDea->starts);
1982
3040
  return MemFree(pCDea);
1983
3041
}
1984
3042
 
1985
 
void           CddExpAlignAlloc(CddExpAlignPtr pCDea, Int4 iLength)
 
3043
void CddExpAlignAlloc(CddExpAlignPtr pCDea, Int4 iLength)
1986
3044
{
1987
3045
  Int4    i;
1988
3046
  
1991
3049
  if (NULL != pCDea->adata) free(pCDea->adata);
1992
3050
  pCDea->adata = MemNew(sizeof(Int4)*pCDea->length);  
1993
3051
  for (i=0;i<pCDea->length;i++) pCDea->adata[i] = -1;
 
3052
  if (NULL != pCDea->starts) free(pCDea->starts);
 
3053
  pCDea->starts = MemNew(sizeof(Int4)*pCDea->length);  
 
3054
  for (i=0;i<pCDea->length;i++) pCDea->starts[i] = 0;
1994
3055
}
1995
3056
 
1996
3057
/*---------------------------------------------------------------------------*/
1997
3058
/* Convert a SeqAlign (pairwise, DenseDiag) into an ExplicitAlignmentObject  */
1998
3059
/*---------------------------------------------------------------------------*/
1999
 
CddExpAlignPtr SeqAlignToCddExpAlign(SeqAlignPtr salp, SeqEntryPtr sep)
 
3060
CddExpAlignPtr LIBCALL SeqAlignToCddExpAlign(SeqAlignPtr salp, SeqEntryPtr sep)
2000
3061
{
2001
3062
  CddExpAlignPtr   pCDea;
2002
3063
  BioseqPtr        bsp1;
2018
3079
      }
2019
3080
      pCDea->adata[ddp->starts[0]+i]=ddp->starts[1]+i;
2020
3081
    }
 
3082
    pCDea->starts[ddp->starts[0]] = 1;
2021
3083
    ddp = ddp->next;
2022
3084
  }
2023
3085
  return(pCDea);
2043
3105
    if (pCDea->adata[i] != -1) {
2044
3106
      pCDeaNew->adata[pCDea->adata[i]] = i;
2045
3107
    }
 
3108
    if (pCDea->starts[i] > 0) {
 
3109
      pCDeaNew->starts[pCDea->adata[i]] = 1;
 
3110
    }
2046
3111
  }
2047
3112
  return(pCDeaNew);
2048
3113
}
2054
3119
{
2055
3120
  SeqAlignPtr   salp;
2056
3121
  DenseDiagPtr  ddp, ddplast = NULL;
2057
 
  Int4          i, last;
 
3122
  Int4          i;
2058
3123
  
2059
3124
  if (!pCDea->ids || !pCDea->ids->next)
2060
3125
   CddSevError("Missing Sequence ID in CddExpAlignToSeqAlign");
2076
3141
      ddp->starts[1] = pCDea->adata[i];
2077
3142
      ddp->len = 1;
2078
3143
      while (i<(pCDea->length-1) && pCDea->adata[i+1]==(pCDea->adata[i]+1) &&
2079
 
             (NULL == iBreakAfter || iBreakAfter[i] == 0)) {
 
3144
             (NULL == iBreakAfter || iBreakAfter[i] == 0) &&
 
3145
             pCDea->starts[i+1] == 0) {
2080
3146
        ddp->len++;
2081
3147
        i++;
2082
3148
      }
2088
3154
      ddplast = ddp;
2089
3155
    }
2090
3156
  }
 
3157
  if (!salp->segs) salp = SeqAlignFree(salp);
2091
3158
  return(salp);
2092
3159
}
2093
3160
 
2121
3188
                iInner, iOuter);
2122
3189
        CddSevError("Error while reindexing explicit alignments!");
2123
3190
      } else  pCDea->adata[iP1] = iP2;
 
3191
      if (pCDea1->starts[i] > 0) pCDea->starts[iP1] = 1;
2124
3192
    }
2125
3193
  }
2126
3194
  if (sip1 && sip2) {
2172
3240
/* return a list of indices for the n most dissimilar sequences in a CD      */
2173
3241
/*---------------------------------------------------------------------------*/
2174
3242
/*---------------------------------------------------------------------------*/
2175
 
Int4Ptr CddMostDiverse(TrianglePtr pTri, Int4 length)
 
3243
Int4Ptr CddMostDiverse(TrianglePtr pTri, Int4 length, Int4 maxdiv)
2176
3244
{
2177
3245
  Int4Ptr     retlist;
2178
3246
  Int4        index, winner, i, j, ncomp;
2179
3247
  FloatHi     sumcomp, mincomp;
2180
3248
  FloatHi     **iMat;
2181
3249
  ScorePtr    psc;
 
3250
  ValNodePtr  vnp;
2182
3251
  
2183
3252
  if (!pTri) return NULL;
2184
3253
  if (length >= pTri->nelements) return NULL;
 
3254
  if (length > maxdiv) length = maxdiv;
 
3255
  retlist = MemNew(length * sizeof(Int4));
 
3256
  retlist[0] = 0;
 
3257
 
 
3258
  if (pTri->div_ranks) {
 
3259
    vnp = pTri->div_ranks; i = 0;
 
3260
    while (vnp) {
 
3261
      retlist[i] = vnp->data.intvalue;
 
3262
      i++; 
 
3263
      if (i >= length) break;
 
3264
      vnp = vnp->next;
 
3265
    }
 
3266
    return(retlist);
 
3267
  }
 
3268
  
2185
3269
  iMat = (FloatHi **) calloc(pTri->nelements,sizeof (FloatHi *));
2186
3270
  for (i=0;i<pTri->nelements;i++) {
2187
3271
    iMat[i] = (FloatHi *) calloc(pTri->nelements,sizeof(FloatHi));
2195
3279
    }
2196
3280
  }
2197
3281
  
2198
 
  retlist = MemNew(length * sizeof(Int4));
2199
 
  retlist[0] = 0;
2200
 
  index = 0; for (index = 1; index < length; index++) {
 
3282
  for (index = 1; index < length; index++) {
2201
3283
    mincomp = 100.0, winner = -1;
2202
3284
    for (i=1;i<pTri->nelements;i++) {
2203
3285
      sumcomp = 0.0; ncomp = 0;
2551
3633
   return(pscHead);
2552
3634
}
2553
3635
 
 
3636
 
 
3637
/*---------------------------------------------------------------------------*/
 
3638
/*---------------------------------------------------------------------------*/
 
3639
/* cloned from salpedit.c and salpacc.c                                      */
 
3640
/*---------------------------------------------------------------------------*/
 
3641
/*---------------------------------------------------------------------------*/
 
3642
SeqAnnotPtr LIBCALL CddSeqAnnotForSeqAlign (SeqAlignPtr salp)
 
3643
{
 
3644
  SeqAnnotPtr sap;
 
3645
 
 
3646
  if (salp==NULL)
 
3647
     return NULL;
 
3648
  sap = SeqAnnotNew ();
 
3649
  sap->type = 2;
 
3650
  sap->data = (Pointer) salp;
 
3651
  return sap;
 
3652
}
 
3653
 
 
3654
SeqAlignPtr LIBCALL CddSeqAlignDup (SeqAlignPtr salp)
 
3655
{
 
3656
  SeqAnnotPtr sap, 
 
3657
              sap2;
 
3658
  SeqAlignPtr salp2 = NULL,
 
3659
              next; 
 
3660
 
 
3661
  next = salp->next;
 
3662
  salp->next = NULL;
 
3663
  sap = CddSeqAnnotForSeqAlign (salp);
 
3664
  sap2 = (SeqAnnotPtr) AsnIoMemCopy ((Pointer) sap, (AsnReadFunc) SeqAnnotAsnRead, (AsnWriteFunc) SeqAnnotAsnWrite);
 
3665
  if (sap2!=NULL) {
 
3666
     salp2 = (SeqAlignPtr)sap2->data;
 
3667
     sap2->data = NULL;
 
3668
     SeqAnnotFree (sap2);
 
3669
  }
 
3670
  if (next != NULL)
 
3671
     salp->next = next;
 
3672
  sap->data = NULL;
 
3673
  SeqAnnotFree (sap);
 
3674
  return salp2;
 
3675
}
 
3676
 
 
3677
 
 
3678
 
 
3679
/*---------------------------------------------------------------------------*/
 
3680
/*---------------------------------------------------------------------------*/
 
3681
/* Duplicate a whole Sequence Alignment (i.e. a linked list of Seqaligns     */
 
3682
/*---------------------------------------------------------------------------*/
 
3683
/*---------------------------------------------------------------------------*/
 
3684
 
 
3685
 
 
3686
SeqAlignPtr LIBCALL SeqAlignSetDup(SeqAlignPtr salp)
 
3687
{
 
3688
  SeqAlignPtr   salpHead = NULL;
 
3689
  SeqAlignPtr   salpTail = NULL;
 
3690
  SeqAlignPtr   salpThis;
 
3691
 
 
3692
  while (salp) {
 
3693
    salpThis = (SeqAlignPtr) CddSeqAlignDup(salp);
 
3694
    if (salpTail) {
 
3695
      salpTail->next = salpThis;
 
3696
      salpTail = salpThis;
 
3697
      salpThis->next = NULL;
 
3698
    } else {
 
3699
      salpHead = salpThis;
 
3700
      salpTail = salpThis;
 
3701
      salpThis->next = NULL;
 
3702
    }
 
3703
    salp = salp->next;
 
3704
  }
 
3705
  return (salpHead);
 
3706
}
 
3707
 
 
3708
 
 
3709
/*---------------------------------------------------------------------------*/
 
3710
/*---------------------------------------------------------------------------*/
 
3711
/* fill internal gaps in a sequence alignment to prepare it for PSSM calc.   */
 
3712
/*---------------------------------------------------------------------------*/
 
3713
/*---------------------------------------------------------------------------*/
 
3714
void LIBCALL CddDegapSeqAlign(SeqAlignPtr salp)
 
3715
{
 
3716
  DenseDiagPtr    ddp;
 
3717
  Int4            iM, iS, iN, iC;
 
3718
  Boolean         bFirst, bInit = TRUE;
 
3719
  BioseqPtr       bsp;
 
3720
  SeqIdPtr        sip;
 
3721
  Int4            mLen, sLen;
 
3722
  
 
3723
  while (salp) {
 
3724
    ddp = salp->segs;
 
3725
    if (bInit) {
 
3726
      sip = ddp->id;
 
3727
      bsp = BioseqLockById(sip);
 
3728
      mLen = bsp->length;
 
3729
      BioseqUnlock(bsp);
 
3730
      bInit = FALSE;
 
3731
    }
 
3732
    sip = ddp->id->next;
 
3733
    bsp = BioseqLockById(sip);
 
3734
    sLen = bsp->length;
 
3735
    BioseqUnlock(bsp);
 
3736
    bFirst = TRUE;
 
3737
    while (ddp) {
 
3738
      if (bFirst) { /* extend first aligned block to N-terminus */
 
3739
        iM = ddp->starts[0];
 
3740
        iS = ddp->starts[1];
 
3741
        if (iM > iS) iM = iS;
 
3742
        ddp->starts[0] -= iM;
 
3743
        ddp->starts[1] -= iM;
 
3744
        ddp->len += iM;
 
3745
        bFirst = FALSE;
 
3746
      }
 
3747
      if (ddp->next) {
 
3748
        iM = ddp->next->starts[0] - (ddp->starts[0]+ddp->len);
 
3749
        iS = ddp->next->starts[1] - (ddp->starts[1]+ddp->len);
 
3750
        if (iM > iS) iM = iS;
 
3751
        if (iM > 0) {
 
3752
          if (iM == 1) {
 
3753
            iN = 1; iC = 0;
 
3754
          } else {
 
3755
            if ((iM % 2) == 0) {
 
3756
              iN = iC = iM / 2;
 
3757
            } else {
 
3758
              iC = iM / 2; iN = iC + 1;
 
3759
            }
 
3760
          }
 
3761
          ASSERT(iN+iC == iM);
 
3762
          ddp->len += iN;
 
3763
          ddp->next->starts[0] -= iC;
 
3764
          ddp->next->starts[1] -= iC;
 
3765
          ddp->next->len += iC;
 
3766
        }
 
3767
      } else { /* extend last aligned block to C-terminus */
 
3768
        iM = mLen - (ddp->starts[0]+ddp->len);      
 
3769
        iS = sLen - (ddp->starts[1]+ddp->len);      
 
3770
        if (iM > iS) iM = iS;
 
3771
        ddp->len += iM;
 
3772
      }
 
3773
      ddp = ddp->next;
 
3774
    }
 
3775
    salp = salp->next;
 
3776
  }
 
3777
}
 
3778
 
 
3779
/*---------------------------------------------------------------------------*/
 
3780
/*---------------------------------------------------------------------------*/
 
3781
/* Duplicate a SeqId, pick the GI if encountering a linked list of several   */
 
3782
/*---------------------------------------------------------------------------*/
 
3783
/*---------------------------------------------------------------------------*/
 
3784
SeqIdPtr LIBCALL CddSeqIdDupPDBGI(SeqIdPtr sipold)
 
3785
{
 
3786
  SeqIdPtr   sip, sipgi = NULL, sippdb = NULL;
 
3787
  
 
3788
  sip = sipold;
 
3789
  while (sip) {
 
3790
    if (sip->choice == SEQID_GI)   sipgi = sip;    
 
3791
    if (sip->choice == SEQID_PDB) sippdb = sip;
 
3792
    sip = sip->next;
 
3793
  }
 
3794
  if (sippdb) return(SeqIdDup(sippdb));
 
3795
  if (sipgi) return(SeqIdDup(sipgi));
 
3796
  return (SeqIdDup(sipold));
 
3797
}
 
3798
 
2554
3799
/*---------------------------------------------------------------------------*/
2555
3800
/*---------------------------------------------------------------------------*/
2556
3801
/* Calculate a weighted 50/50 consensus sequence and make it new master      */
2572
3817
  DenseDiagPtr        ddp;
2573
3818
  Int4Ptr             trunc_on_virtual;
2574
3819
  Int4Ptr             trunc_aligned;
2575
 
  Int4Ptr             virtual_on_trunc;
2576
3820
  ObjectIdPtr         oidp;
2577
3821
  SeqAlignPtr         salpThis, salpNew, salpTail;
2578
3822
  SeqEntryPtr         sepNew;
2579
 
  SeqIdPtr            sipThis;
 
3823
  SeqIdPtr            sipThis, sipCopy;
2580
3824
  SeqMapTablePtr      smtp;
2581
3825
  SeqPortPtr          spp;
2582
3826
  Uint1Ptr            buffer;
2586
3830
  Int4                offset, i, j , k, ipt, seqpos;
2587
3831
  Int4                m_last, m_frst, s_insert, m_end;
2588
3832
  Int4                s_last, s_frst, maxpos = 0;
2589
 
  Uint1               aatype;
2590
3833
  
2591
3834
  if (!salp) return (NULL);
2592
3835
  if (salp->type != SAT_PARTIAL || salp->segtype != SAS_DENDIAG) return(NULL);
2702
3945
  i = 0;
2703
3946
  salpThis = salp;
2704
3947
  ddp = salpThis->segs;
2705
 
  sipThis = ddp->id;
2706
 
  bspThis = CddRetrieveBioseqById(sipThis,sep);
 
3948
  sipCopy = SeqIdDup(ddp->id); sipCopy->next = NULL;
 
3949
  bspThis = CddRetrieveBioseqById(sipCopy,sep);
 
3950
  SeqIdFree(sipCopy);
2707
3951
  pCAW[i].bsp = bspThis;
2708
3952
  spp = SeqPortNew(bspThis, 0, LAST_RESIDUE, Seq_strand_unknown,
2709
3953
                   Seq_code_ncbistdaa);
2717
3961
    ASSERT (k>=0 && k<maxextlen);
2718
3962
    pCEACell[i][k].seqpos = seqpos;
2719
3963
    pCEACell[i][k].aatype = buffer[seqpos];
2720
 
    if (pCEACell[i][k].aatype < 0 || pCEACell[i][k].aatype > 26) {
 
3964
    if (pCEACell[i][k].aatype > (Uint1) 26) {
2721
3965
      pCEACell[i][k].aatype = 26;
2722
3966
    }
2723
3967
  }
2744
3988
        ASSERT (k>=0 && k<maxextlen);
2745
3989
        pCEACell[i][k].seqpos = ddp->starts[1] + j - ddp->starts[0];
2746
3990
        pCEACell[i][k].aatype = buffer[pCEACell[i][k].seqpos];
2747
 
        if (pCEACell[i][k].aatype < 0 || pCEACell[i][k].aatype > 26) {
 
3991
        if (pCEACell[i][k].aatype > (Uint1) 26) {
2748
3992
          pCEACell[i][k].aatype = 26;
2749
3993
        }
2750
3994
      }
2760
4004
            ASSERT (k>=0 && k<maxextlen);
2761
4005
            pCEACell[i][k].seqpos = j;
2762
4006
            pCEACell[i][k].aatype = buffer[j];
2763
 
            if (pCEACell[i][k].aatype < 0 || pCEACell[i][k].aatype > 26) {
 
4007
            if (pCEACell[i][k].aatype > (Uint1) 26) {
2764
4008
              pCEACell[i][k].aatype = 26;
2765
4009
            }
2766
4010
          }
2769
4013
            ASSERT (k>=0 && k<maxextlen);
2770
4014
            pCEACell[i][k].seqpos = j;
2771
4015
            pCEACell[i][k].aatype = buffer[j];
2772
 
            if (pCEACell[i][k].aatype < 0 || pCEACell[i][k].aatype > 26) {
 
4016
            if (pCEACell[i][k].aatype > (Uint1) 26) {
2773
4017
              pCEACell[i][k].aatype = 26;
2774
4018
            }
2775
4019
          }     
2816
4060
    if (pCEACol[k].occpos >= maxpos) {
2817
4061
      ContWeight++;
2818
4062
      for (i=0;i<nalign;i++) {
2819
 
        if (pCEACell[i][k].seqpos != -1 && pCEACell[i][k].aatype < 26) {
 
4063
        if (pCEACell[i][k].seqpos != -1 && pCEACell[i][k].aatype < (Uint1) 26) {
2820
4064
          pCAW[i].nposaligned++;
2821
4065
          if (pCEACol[k].typecount[pCEACell[i][k].aatype] > 0) {
2822
4066
            pCAW[i].weight += 1.0 / (pCEACol[k].ntypes * 
2845
4089
  ConsensusLen = 0;
2846
4090
  for (k=0;k<maxextlen;k++) {
2847
4091
    for (i=0;i<nalign;i++) {
2848
 
      if (pCEACell[i][k].seqpos != -1 && pCEACell[i][k].aatype < 26) {
 
4092
      if (pCEACell[i][k].seqpos != -1 && pCEACell[i][k].aatype < (Uint1) 26) {
2849
4093
        pCEACol[k].w_occpos += pCAW[i].weight;
2850
4094
        pCEACol[k].wtypefreq[pCEACell[i][k].aatype]+=pCAW[i].weight;
2851
4095
      }
2894
4138
  bspThis = pCAW[0].bsp;
2895
4139
  pCDea = CddExpAlignNew();
2896
4140
  CddExpAlignAlloc(pCDea,bspThis->length);
2897
 
  pCDea->ids = SeqIdDup(bspThis->id);
 
4141
  pCDea->ids = CddSeqIdDupPDBGI(bspThis->id);
2898
4142
  pCDea->ids->next = SeqIdDup((*bspCons)->id);
2899
4143
  pCDea->bIdAlloc = TRUE;
2900
4144
  for (k=0;k<maxextlen;k++) {
2902
4146
      pCDea->adata[pCEACell[0][k].seqpos]=pCEACol[k].conpos;
2903
4147
    }
2904
4148
  }
 
4149
/*---------------------------------------------------------------------------*/
 
4150
/* record block structure present in the input alignment to prevent merges   */
 
4151
/*---------------------------------------------------------------------------*/
 
4152
  ddp = salp->segs;
 
4153
  while (ddp) {
 
4154
    pCDea->starts[ddp->starts[0]] = 1;
 
4155
    ddp = ddp->next;
 
4156
  }
2905
4157
  *salpCons = CddExpAlignToSeqAlign(pCDea,NULL);
2906
4158
  pCDea = CddExpAlignFree(pCDea);
2907
4159
 
2922
4174
    }
2923
4175
    salpThis = CddExpAlignToSeqAlign(pCDea,NULL);
2924
4176
    pCDea = CddExpAlignFree(pCDea);
2925
 
    if (!salpNew) {
2926
 
      salpNew = salpThis;
2927
 
    } else {
2928
 
      salpTail->next = salpThis;
 
4177
    if (salpThis) {
 
4178
      if (!salpNew) {
 
4179
        salpNew = salpThis;
 
4180
      } else {
 
4181
        salpTail->next = salpThis;
 
4182
      }
 
4183
      salpTail = salpThis;
2929
4184
    }
2930
 
    salpTail = salpThis;
2931
4185
  }
2932
4186
 
2933
4187
/*---------------------------------------------------------------------------*/
2949
4203
}
2950
4204
 
2951
4205
 
2952
 
static void CddCposCompPart1(SeqAlignPtr listOfSeqAligns, BioseqPtr bsp,
2953
 
                             compactSearchItems* compactSearch, ValNodePtr* LetterHead,
2954
 
                             posSearchItems* posSearch) {
 
4206
/*---------------------------------------------------------------------------*/
2955
4207
/*---------------------------------------------------------------------------*/
2956
4208
/* first part of CddCposComp()                                               */
2957
4209
/* basically returns compactSearch, LetterHead, and posSearch                */
2958
4210
/*---------------------------------------------------------------------------*/
 
4211
/*---------------------------------------------------------------------------*/
 
4212
static void CddCposCompPart1(SeqAlignPtr listOfSeqAligns, BioseqPtr bsp,
 
4213
                             compactSearchItems* compactSearch, ValNodePtr* LetterHead,
 
4214
                             posSearchItems* posSearch) {
2959
4215
    Int4                numalign, numseq;    /*number of alignments and seqs */
2960
4216
    BLAST_ScoreBlkPtr   sbp;
2961
4217
    SeqLocPtr           private_slp = NULL;
3268
4524
{
3269
4525
    PGPBlastOptionsPtr  bop;
3270
4526
    BLAST_OptionsBlkPtr options;
3271
 
    SeqEntryPtr         sep;
3272
 
    Boolean             is_dna;
3273
 
    ObjectIdPtr         obidp;
3274
4527
 
3275
4528
    bop = MemNew(sizeof(PGPBlastOptions));
3276
4529
    bop->blast_database   = StringSave("nr");
3416
4669
  return(TRUE);
3417
4670
}
3418
4671
 
 
4672
/*---------------------------------------------------------------------------*/
 
4673
/*---------------------------------------------------------------------------*/
 
4674
/*---------------------------------------------------------------------------*/
 
4675
/*---------------------------------------------------------------------------*/
 
4676
static void BlastKarlinBlkCopy(BLAST_KarlinBlkPtr kbp_in,BLAST_KarlinBlkPtr kbp_out)
 
4677
{ 
 
4678
  kbp_out->Lambda = kbp_in->Lambda;
 
4679
  kbp_out->K = kbp_in->K;
 
4680
  kbp_out->logK = kbp_in->logK;
 
4681
  kbp_out->H = kbp_in->H;
 
4682
/*  kbp_out->Lambda_real = kbp_in->Lambda_real;
 
4683
  kbp_out->K_real = kbp_in->K_real;
 
4684
  kbp_out->logK_real = kbp_in->logK_real;
 
4685
  kbp_out->H_real = kbp_in->H_real;
 
4686
  kbp_out->q_frame = kbp_in->q_frame;
 
4687
  kbp_out->s_frame = kbp_in->s_frame; */
 
4688
  kbp_out->paramC = kbp_in->paramC;
 
4689
}
3419
4690
 
3420
4691
/*---------------------------------------------------------------------------*/
3421
4692
/*---------------------------------------------------------------------------*/
3422
 
/* updated version of the PSSM calculation routine, modified from blastpgp.c */
3423
 
/* as of 12/27/2000                                                          */
 
4693
/* version of the PSSM calculation routine which does not return Karlin-     */
 
4694
/* Altschul parameters                                                       */
3424
4695
/*---------------------------------------------------------------------------*/
3425
4696
/*---------------------------------------------------------------------------*/
3426
4697
Seq_Mtf * LIBCALL CddDenDiagCposComp2(BioseqPtr bspFake, Int4 iPseudo,
3428
4699
                                      BioseqPtr bspOut, double Weight,
3429
4700
                                      double ScaleFactor, CharPtr matrix_name)
3430
4701
{
 
4702
  return CddDenDiagCposComp2KBP(bspFake, iPseudo, salp, pcdd, bspOut,
 
4703
                                Weight, ScaleFactor, matrix_name, NULL);
 
4704
}
 
4705
 
 
4706
/*---------------------------------------------------------------------------*/
 
4707
/*---------------------------------------------------------------------------*/
 
4708
/* updated version of the PSSM calculation routine, modified from blastpgp.c */
 
4709
/* as of 12/27/2000                                                          */
 
4710
/*---------------------------------------------------------------------------*/
 
4711
/*---------------------------------------------------------------------------*/
 
4712
Seq_Mtf * LIBCALL CddDenDiagCposComp2KBP(BioseqPtr bspFake, Int4 iPseudo,
 
4713
                                         SeqAlignPtr salp, CddPtr pcdd,
 
4714
                                         BioseqPtr bspOut, double Weight,
 
4715
                                         double ScaleFactor, CharPtr matrix_name,
 
4716
                                         BLAST_KarlinBlkPtr kbp)
 
4717
{
3431
4718
  PGPBlastOptionsPtr      bop;
3432
4719
  BlastSearchBlkPtr       search;
3433
4720
  posSearchItems         *posSearch = NULL;
3435
4722
  ValNodePtr              ColumnHead, LetterHead, newRow, RowHead, newLetter;
3436
4723
  SeqCodeTablePtr         sctp;
3437
4724
  Int4                    numSeqs, numASeqs, alignLength, index, a, iNew, i;
3438
 
  Int4                    scale_factor, j, OutLen, jj, c;
 
4725
  Int4                    scale_factor, j, OutLen, jj;
3439
4726
  Char                    ckptFileName[PATH_MAX];
3440
4727
  Char                    cseqFileName[PATH_MAX];
3441
4728
  Char                    mtrxFileName[PATH_MAX];
3442
4729
  FILE                   *fp;
3443
4730
  Boolean                 bHasConsensus;
 
4731
  Boolean                 bWriteOut = FALSE;
3444
4732
  Nlm_FloatHi            *AInf, SumAInf = 0.0;
3445
4733
  char*                   Input;    /* order of residue-types from CD's      */
3446
4734
  char*                   Output;   /*order of res-types needed for threading*/
3447
4735
  Boolean*                Coverage; /* for making sure all columns are filled*/
3448
4736
  Seq_Mtf*                pssm;
3449
4737
  CharPtr                 cAccession;
 
4738
  Boolean                 scalingOkay = TRUE;
3450
4739
 
3451
4740
  cAccession = CddGetAccession(pcdd);
3452
4741
  if (pcdd) {
3453
4742
    bHasConsensus = CddHasConsensus(pcdd);
 
4743
    bWriteOut = CddHasDescription(pcdd,"Write out matrix");
3454
4744
  } else {                                        /* determine from seqalign */
3455
4745
    bHasConsensus = SeqAlignHasConsensus(salp);
3456
4746
  }
3457
4747
  if (iPseudo <= 0) {                             /* need to determine first */
3458
 
    AInf = SeqAlignInform(salp, bspFake, bHasConsensus); 
 
4748
    AInf = SeqAlignInform(salp, bspFake, bHasConsensus, 0); 
3459
4749
    for (i=0;i<bspFake->length;i++) SumAInf += AInf[i];
3460
 
    if (SumAInf > 84) iPseudo = 10;                    /* purely empirical   */
3461
 
    else if (SumAInf > 55) iPseudo = 7;
3462
 
    else if (SumAInf > 43) iPseudo = 5;
3463
 
    else if (SumAInf > 41.5) iPseudo = 4;
3464
 
    else if (SumAInf > 40) iPseudo = 3;
3465
 
    else if (SumAInf > 39) iPseudo = 2;
3466
 
    else iPseudo = 1;
 
4750
    if      (SumAInf > 84  ) iPseudo = 10;             /* purely empirical   */
 
4751
    else if (SumAInf > 55  ) iPseudo =  7;
 
4752
    else if (SumAInf > 43  ) iPseudo =  5;
 
4753
    else if (SumAInf > 41.5) iPseudo =  4;
 
4754
    else if (SumAInf > 40  ) iPseudo =  3;
 
4755
    else if (SumAInf > 39  ) iPseudo =  2;
 
4756
    else                     iPseudo =  1;
3467
4757
    MemFree(AInf);
3468
 
    if (pcdd) printf("%s AInf:%6.1f PseudoCt: %d\n",cAccession,SumAInf, iPseudo);
 
4758
    if (pcdd && bWriteOut) printf("%s AInf:%6.1f PseudoCt: %d\n",cAccession,SumAInf, iPseudo);
3469
4759
  }
3470
4760
  if (NULL == matrix_name) {
3471
4761
    static CharPtr defaultMatrix = "BLOSUM62";
3481
4771
    search = BLASTSetUpSearch(bop->fake_bsp, "blastp",
3482
4772
                              bop->query_bsp->length,
3483
4773
                              0, NULL, bop->options, NULL);
3484
 
  
3485
4774
  }
3486
 
 
3487
4775
  posSearch = NULL;
3488
4776
  compactSearch = NULL;
3489
4777
  posSearch = (posSearchItems *)MemNew(sizeof(posSearchItems));
3490
 
 
3491
4778
  compactSearch = compactSearchNew(compactSearch);
3492
4779
  copySearchItems(compactSearch, search, bop->options->matrix);
3493
4780
  posInitializeInformation(posSearch,search);
3517
4804
/*---------------------------------------------------------------------------*/
3518
4805
/* Construct name for checkpoint file and write out (if in a CD context)     */
3519
4806
/*---------------------------------------------------------------------------*/
3520
 
  if (pcdd) {
 
4807
  if (pcdd && bWriteOut) {
3521
4808
    strcpy(ckptFileName,cAccession);
3522
4809
    strcat(ckptFileName,CKPTEXT);
3523
4810
    if (NULL != ckptFileName)
3524
4811
      posTakeCheckpoint(posSearch, compactSearch, ckptFileName, NULL);
3525
4812
  }
3526
4813
 
3527
 
  posFreqsToMatrix(posSearch,compactSearch, NULL, 1);
 
4814
  posFreqsToMatrix(posSearch,compactSearch);
3528
4815
/*---------------------------------------------------------------------------*/
3529
4816
/* kludge to set the matrix column for matches to * (stop codons) to the     */
3530
4817
/* most negative value used, so that matrix-file corresponds to what makemat */
3531
4818
/* produces.                                                                 */
3532
4819
/*---------------------------------------------------------------------------*/
3533
4820
  ASSERT(compactSearch->alphabetSize == 26);                 /* sanity check */
 
4821
/*
3534
4822
  for (c=0;c<compactSearch->qlength;c++) {
3535
4823
    posSearch->posPrivateMatrix[c][compactSearch->alphabetSize-1] = BLAST_SCORE_MIN;
3536
4824
  }
3537
 
 
3538
 
  posScaling(posSearch, compactSearch);
 
4825
*/
 
4826
/*  posScaling(posSearch, compactSearch); */
 
4827
  scalingOkay = impalaScaling(posSearch, compactSearch, 1.0, TRUE);
 
4828
  if (!scalingOkay) goto cleanup;
3539
4829
 
3540
4830
  if (pcdd) {
3541
4831
/*---------------------------------------------------------------------------*/
3542
4832
/* also, write out the PSSM as an ASCII-formatted MATRIX file, for use with  */
3543
4833
/* copymat                                                                   */
3544
4834
/*---------------------------------------------------------------------------*/
3545
 
    strcpy(mtrxFileName,cAccession);
 
4835
/*    strcpy(mtrxFileName,cAccession);
3546
4836
    strcat(mtrxFileName,MTRXEXT);
3547
4837
    CddtakeMatrixCheckpoint(compactSearch,posSearch,search->sbp,mtrxFileName,
3548
4838
                            NULL,FALSE,1.0);
3549
 
 
 
4839
*/
 
4840
/*---------------------------------------------------------------------------*/
 
4841
/* in compliance with makemat, set the score for the X-column                */
 
4842
/*---------------------------------------------------------------------------*/
 
4843
    for (index = 0; index<compactSearch->qlength;index++) {
 
4844
      posSearch->posMatrix[index][21] = Xscore;
 
4845
    }
3550
4846
    sctp = SeqCodeTableFind(Seq_code_ncbistdaa);
3551
4847
    LetterHead = NULL;
3552
4848
    for (a=0;a<compactSearch->alphabetSize;a++) {
3617
4913
    MemFree(Coverage);
3618
4914
  }
3619
4915
 
3620
 
  if (bspOut && pcdd) {
 
4916
  if (bspOut && pcdd && bWriteOut) {
3621
4917
    strcpy(cseqFileName,cAccession);
3622
4918
    strcat(cseqFileName,CSEQEXT);
3623
4919
    fp = FileOpen(cseqFileName, "w");
3631
4927
/*---------------------------------------------------------------------------*/
3632
4928
/* clean up                                                                  */
3633
4929
/*---------------------------------------------------------------------------*/
 
4930
cleanup:
3634
4931
  MemFree(cAccession);
3635
4932
  posSearch->NumSequences = posSearch->posNumSequences;
3636
4933
  posSearch->QuerySize = alignLength;
3637
4934
  CddposFreeMemory(posSearch);
3638
4935
  MemFree(posSearch);
 
4936
  if (kbp) {
 
4937
    BlastKarlinBlkCopy(*(compactSearch->kbp_gap_psi),kbp);
 
4938
  }
3639
4939
  CddposFreeMemory2(compactSearch);
3640
4940
  MemFree(compactSearch);
3641
4941
  search = BlastSearchBlkDestruct(search);
3642
4942
  bop->options = BLASTOptionDelete(bop->options);
3643
4943
  MemFree(bop->blast_database);
3644
4944
  MemFree(bop);
3645
 
  if (NULL == pcdd) {        /* return pssm if used from within the threader */
 
4945
  if (NULL == pcdd && scalingOkay) {  /* return pssm if from within threader */
3646
4946
    return(pssm);
3647
4947
  }
 
4948
  if (NULL != pcdd && !scalingOkay) { /* fields are empty if scaling failed  */
 
4949
    pcdd->posfreq = NULL;
 
4950
    pcdd->scoremat = NULL;
 
4951
  }
3648
4952
  return NULL;
3649
4953
}
3650
4954
 
3757
5061
}
3758
5062
 
3759
5063
/*---------------------------------------------------------------------------*/
3760
 
/*---------------------------------------------------------------------------*/
 
5064
/* Translate a SeqLoc into an Integer list of residue numbers                */
 
5065
/*---------------------------------------------------------------------------*/
 
5066
Int4Ptr LIBCALL CddGetFeatLocList(SeqLocPtr location, Int4 *nres)
 
5067
{
 
5068
  SeqLocPtr     slp;
 
5069
  SeqIntPtr     sintp;
 
5070
  PackSeqPntPtr pspp;
 
5071
  SeqPntPtr     spp;
 
5072
  Int4Ptr       presnum = NULL;
 
5073
  Int4          i;
 
5074
  Int4          nSubLoc;
 
5075
  Int4Ptr       pSubLoc;
 
5076
 
 
5077
  presnum = MemNew(100 * sizeof(Int4));
 
5078
  *nres = 0;
 
5079
  if (!location) return NULL;
 
5080
  switch (location->choice) {
 
5081
    case SEQLOC_NULL:
 
5082
    case SEQLOC_FEAT:
 
5083
    case SEQLOC_BOND:
 
5084
    case SEQLOC_EMPTY:
 
5085
    case SEQLOC_WHOLE:
 
5086
      break;
 
5087
    case SEQLOC_INT:
 
5088
      sintp = (SeqIntPtr) location->data.ptrvalue;
 
5089
      for (i=sintp->from;i<=sintp->to;i++) {
 
5090
        presnum[*nres] = i;
 
5091
        (*nres)++; if (*nres >= 100) presnum = Nlm_MemMore(presnum,((*nres)+1)*sizeof(Int4));
 
5092
      }
 
5093
      break;
 
5094
    case SEQLOC_PNT:
 
5095
      spp = location->data.ptrvalue;
 
5096
      presnum[*nres] = spp->point;
 
5097
        (*nres)++; if (*nres >= 100) presnum = Nlm_MemMore(presnum,((*nres)+1)*sizeof(Int4));
 
5098
      break;
 
5099
    case SEQLOC_PACKED_PNT:
 
5100
      pspp = location->data.ptrvalue;
 
5101
      while (pspp) {
 
5102
        for (i=0;i<pspp->used;i++) {
 
5103
          presnum[*nres] = pspp->pnts[i];
 
5104
          (*nres)++; if (*nres >= 100) presnum = Nlm_MemMore(presnum,((*nres)+1)*sizeof(Int4));
 
5105
        }
 
5106
        pspp = pspp->next;
 
5107
      }
 
5108
      break;
 
5109
    case SEQLOC_PACKED_INT:
 
5110
    case SEQLOC_MIX:
 
5111
    case SEQLOC_EQUIV:
 
5112
      slp = (SeqLocPtr)location->data.ptrvalue;
 
5113
      while (slp) {
 
5114
        pSubLoc = CddGetFeatLocList(slp, &nSubLoc);
 
5115
        for (i=0;i<nSubLoc;i++) {
 
5116
          presnum[*nres] = pSubLoc[i];
 
5117
          (*nres)++; if (*nres >= 100) presnum = Nlm_MemMore(presnum,((*nres)+1)*sizeof(Int4));
 
5118
        }       
 
5119
        slp = slp->next;
 
5120
      }      
 
5121
    default:
 
5122
      break;
 
5123
  }
 
5124
  return(presnum);
 
5125
}
 
5126
 
 
5127
/*---------------------------------------------------------------------------*/
 
5128
/* check for violations of an interval defined on the aligned master with    */
 
5129
/* respect to aligned residues. Returns -1 if everything is OK, returns      */
 
5130
/* the first location on the master (counting from 0) if a violation is      */
 
5131
/* detected, returns INT4_MAX if negative interval numbers are encountered   */
 
5132
/*---------------------------------------------------------------------------*/
 
5133
static Int4 CddSeqLocInExpAlign(SeqLocPtr location, CddExpAlignPtr eap)
 
5134
{
 
5135
  SeqIntPtr     sintp;
 
5136
  SeqPntPtr     spp;
 
5137
  PackSeqPntPtr pspp;
 
5138
  SeqLocPtr     slp;
 
5139
  Int4          i;
 
5140
 
 
5141
  if (!location) return(-1);
 
5142
  switch (location->choice) {
 
5143
    case SEQLOC_NULL:
 
5144
    case SEQLOC_FEAT:
 
5145
    case SEQLOC_BOND:
 
5146
    case SEQLOC_EMPTY:
 
5147
    case SEQLOC_WHOLE:
 
5148
      break;
 
5149
    case SEQLOC_INT:
 
5150
      sintp = (SeqIntPtr) location->data.ptrvalue;
 
5151
      if (sintp->from < 0 || sintp->to < 0) return(INT4_MAX);
 
5152
      for (i=sintp->from;i<=sintp->to;i++) {
 
5153
        if (eap->adata[i] < 0) return(i);
 
5154
      }
 
5155
      break;
 
5156
    case SEQLOC_PNT:
 
5157
      spp = location->data.ptrvalue;
 
5158
      if (spp->point < 0) return(INT4_MAX);
 
5159
      if (eap->adata[spp->point] < 0) return(spp->point);
 
5160
      break;
 
5161
    case SEQLOC_PACKED_PNT:
 
5162
      pspp = location->data.ptrvalue;
 
5163
      while (pspp) {
 
5164
        for (i=0;i<pspp->used;i++) {
 
5165
          if (pspp->pnts[i] < 0) return (INT4_MAX);
 
5166
          if (eap->adata[pspp->pnts[i]] < 0) return(pspp->pnts[i]);
 
5167
        }
 
5168
        pspp = pspp->next;
 
5169
      }
 
5170
      break;
 
5171
    case SEQLOC_PACKED_INT:
 
5172
    case SEQLOC_MIX:
 
5173
    case SEQLOC_EQUIV:
 
5174
      slp = (SeqLocPtr)location->data.ptrvalue;
 
5175
      while (slp) {
 
5176
        i = CddSeqLocInExpAlign(slp,eap);
 
5177
        if (i >= 0) return(i);
 
5178
        slp = slp->next;
 
5179
      }      
 
5180
    default:
 
5181
      break;
 
5182
  }
 
5183
  return(-1);
 
5184
}
 
5185
 
 
5186
 
3761
5187
/*---------------------------------------------------------------------------*/
3762
5188
/*---------------------------------------------------------------------------*/
3763
5189
static void CddRelocateSeqLoc(SeqLocPtr location, SeqIdPtr sip, Int4 *ali)
3778
5204
    case SEQLOC_WHOLE:
3779
5205
      SeqIdFree(location->data.ptrvalue);
3780
5206
      location->data.ptrvalue = (SeqIdPtr) SeqIdDup(sip);
 
5207
      break;
3781
5208
    case SEQLOC_INT:
3782
5209
      sintp = (SeqIntPtr) location->data.ptrvalue;
3783
5210
      SeqIdFree(sintp->id);
3784
5211
      sintp->id = SeqIdDup(sip);
3785
5212
      sintp->from = ali[sintp->from];
3786
5213
      sintp->to   = ali[sintp->to];
 
5214
      if (sintp->from < 0 || sintp->to < 0) CddSevError("Invalid SeqLoc Interval!");
3787
5215
      break;
3788
5216
    case SEQLOC_PNT:
3789
5217
      spp = location->data.ptrvalue;
3790
5218
      SeqIdFree(spp->id);
3791
5219
      spp->id = (SeqIdPtr) SeqIdDup(sip);
3792
5220
      spp->point = ali[spp->point];
 
5221
      if (spp->point < 0) CddSevError("Invalid SeqLoc Interval!");
3793
5222
      break;
3794
5223
    case SEQLOC_PACKED_PNT:
3795
5224
      pspp = location->data.ptrvalue;
3798
5227
        pspp->id = (SeqIdPtr) SeqIdDup(sip);
3799
5228
        for (i=0;i<pspp->used;i++) {
3800
5229
          pspp->pnts[i] = ali[pspp->pnts[i]];   
 
5230
          if (pspp->pnts[i] < 0) CddSevError("Invalid SeqLoc Interval!");
3801
5231
        }
3802
5232
        pspp = pspp->next;
3803
5233
      }
3836
5266
{
3837
5267
  CddExpAlignPtr    pCDeaO= NULL, pCDeaN = NULL, pCDeaON = NULL;
3838
5268
  SeqLocPtr         slp;
3839
 
  SeqIdPtr          oldMaster, sip;
 
5269
  SeqIdPtr          oldMaster, sip, sipTmp;
3840
5270
  Boolean           bOldIsTop = FALSE;
3841
5271
  Boolean           bNewIsTop = FALSE;
3842
5272
  DenseDiagPtr      ddp;
3843
5273
  SeqAlignPtr       salpThis;
3844
5274
  BioseqPtr         bsp;
 
5275
  SeqEntryPtr       sep;
3845
5276
  AlignAnnotPtr     aap;
3846
5277
  
3847
5278
  slp = (SeqLocPtr) oldannot->location;
3848
5279
  oldMaster = CddFindSeqIdInAlignAnnot(oldannot);
 
5280
  if (oldMaster->choice == SEQID_GI) {
 
5281
    sep = bssp->seq_set; while (sep) {
 
5282
      if (sep->choice == 1) {
 
5283
        bsp = sep->data.ptrvalue;
 
5284
        if (CddSameSip(oldMaster, bsp->id)) {
 
5285
          sipTmp = bsp->id;
 
5286
          while (sipTmp) {
 
5287
            if (sipTmp->choice == SEQID_PDB) {
 
5288
              oldMaster = sipTmp;
 
5289
              break;      
 
5290
            }
 
5291
            sipTmp = sipTmp->next;
 
5292
          }
 
5293
          break;
 
5294
        }
 
5295
      }
 
5296
      sep = sep->next;
 
5297
    }
 
5298
  }
3849
5299
  ddp = salp->segs;
3850
5300
  sip = SeqIdDup(ddp->id);  SeqIdFree(sip->next);
3851
5301
  if (CddSameSip(oldMaster,sip)) bOldIsTop = TRUE;
3883
5333
  if (bNewIsTop) {                             /* need to revert alignment  */
3884
5334
    pCDeaON = InvertCddExpAlign(pCDeaON,bssp->seq_set);
3885
5335
  }
 
5336
  if (!pCDeaON) CddSevError("Error in CDDTransferAlignAnnot, can not locate alignment data!");
3886
5337
  aap = oldannot; while (aap) {
3887
5338
    slp = aap->location;
3888
5339
    while (slp) {
3949
5400
/********************************************************************/
3950
5401
void LIBCALL CddShrinkBioseq(BioseqPtr bsp)
3951
5402
{
3952
 
    SeqDescrPtr next,sdp,sdp2=NULL;
3953
 
    SeqAnnotPtr sap,sap_next,sap2=NULL;
3954
 
 
3955
 
    sdp=bsp->descr;
3956
 
    while (sdp != NULL) {
3957
 
        next = sdp->next;
3958
 
        sdp->next = NULL;
3959
 
        if(sdp->choice == Seq_descr_title ||
3960
 
           sdp->choice == Seq_descr_org ||
3961
 
           sdp->choice == Seq_descr_source) {
3962
 
            ValNodeLink(&sdp2, sdp);
3963
 
        } else SeqDescFree(sdp);
3964
 
        sdp = next;
3965
 
    }
3966
 
    bsp->descr = sdp2;
3967
 
 
3968
 
    sap = bsp->annot;
3969
 
    while (sap != NULL) {
3970
 
        sap_next = sap->next;
3971
 
        if(sap->type == 4 && ((SeqIdPtr)sap->data)->choice == SEQID_GENERAL
3972
 
           && strcmp(((DbtagPtr)((SeqIdPtr)sap->data)->
3973
 
                      data.ptrvalue)->db,"mmdb") == 0 && sap2 == NULL) {
3974
 
            sap2 = sap;
3975
 
        } else SeqAnnotFree(sap);
3976
 
        sap = sap_next;
3977
 
    }
3978
 
    bsp->annot = sap2;
3979
 
    
 
5403
  SeqDescrPtr next, sdp, sdp2              = NULL;
 
5404
  SeqAnnotPtr sap, sap_next, thissap, sap2 = NULL;
 
5405
 
 
5406
  sdp=bsp->descr;
 
5407
  while (sdp != NULL) {
 
5408
    next = sdp->next;
 
5409
    sdp->next = NULL;
 
5410
    if(sdp->choice == Seq_descr_title ||
 
5411
      sdp->choice == Seq_descr_org ||
 
5412
      sdp->choice == Seq_descr_pdb ||
 
5413
      sdp->choice == Seq_descr_source) {
 
5414
        ValNodeLink(&sdp2, sdp);
 
5415
      } else SeqDescFree(sdp);
 
5416
    sdp = next;
 
5417
  }
 
5418
  bsp->descr = sdp2;
 
5419
 
 
5420
  sap = bsp->annot;
 
5421
  while (sap != NULL) {
 
5422
    sap_next = sap->next;
 
5423
    if ((sap->type == 4 && ((SeqIdPtr)sap->data)->choice == SEQID_GENERAL
 
5424
         && Nlm_StrCmp(((DbtagPtr)((SeqIdPtr)sap->data)->data.ptrvalue)->db,"mmdb") == 0)
 
5425
        || (sap->type == 1 && ((SeqFeatPtr)sap->data)->data.choice == SEQFEAT_PROT)) {
 
5426
      if (!sap2) {
 
5427
        sap2 = sap;
 
5428
      } else {
 
5429
        thissap = sap2; while (thissap->next) thissap = thissap->next;
 
5430
        thissap->next = sap;
 
5431
      }
 
5432
      sap->next = NULL;
 
5433
    } else SeqAnnotFree(sap);
 
5434
    sap = sap_next;
 
5435
  }
 
5436
  bsp->annot = sap2;
3980
5437
}
3981
 
/********************************************************************/
3982
 
/* Make a list of m-s seq-aligns proper, aligning only columns      */
3983
 
/* found in all pairs.                                              */
3984
 
/********************************************************************/
 
5438
 
 
5439
/*---------------------------------------------------------------------------*/
 
5440
/*---------------------------------------------------------------------------*/
 
5441
/* Make a list of m-s seq-aligns proper, aligning only columns               */
 
5442
/* found in all pairs.                                                       */
 
5443
/*---------------------------------------------------------------------------*/
 
5444
/*---------------------------------------------------------------------------*/
3985
5445
Int2 LIBCALL CddTrimSeqAligns(CddPtr pcdd)
3986
5446
{
3987
5447
    CddExpAlignPtr      ceap,tceap;
3990
5450
    ValNodePtr          vnp,vhead=NULL;
3991
5451
    Int2                i;
3992
5452
    Int4Ptr             iBreakAfter;
 
5453
    DenseDiagPtr        ddp;
3993
5454
 
3994
5455
    if(pcdd == NULL) return 1;
3995
5456
    sep = ((BioseqSetPtr)pcdd->sequences->data.ptrvalue)->seq_set;
3996
5457
    head = pcdd->seqannot->data;
3997
 
    if(head == NULL || head->next == NULL) return 1;
 
5458
    if(head == NULL /*|| head->next == NULL */) return 1;
3998
5459
 
3999
5460
    tceap = (CddExpAlignPtr)SeqAlignToCddExpAlign(head, sep);
4000
5461
    iBreakAfter = MemNew(tceap->length * sizeof(Int4));
 
5462
/*---------------------------------------------------------------------------*/
 
5463
/* record existing block structure to make sure adjacent blocks aren't merged*/
 
5464
/*---------------------------------------------------------------------------*/
4001
5465
    for (i=0;i<tceap->length;i++) iBreakAfter[i] = 0;
 
5466
    ddp = head->segs;
 
5467
    while (ddp) {
 
5468
      iBreakAfter[ddp->starts[0]+ddp->len-1] = 1;
 
5469
      ddp = ddp->next;
 
5470
    }    
4002
5471
    for(salp=head; salp; salp=salp->next) {
4003
5472
      ceap = (CddExpAlignPtr)SeqAlignToCddExpAlign(salp, sep);
4004
5473
      ValNodeAddPointer(&vhead, 0, ceap);
4133
5602
 
4134
5603
    return 0;
4135
5604
}
 
5605
 
 
5606
/*---------------------------------------------------------------------------*/
 
5607
/* fix character string for export into XML                                  */
 
5608
/*---------------------------------------------------------------------------*/
 
5609
static CharPtr CddFixForXML(CharPtr pc)
 
5610
{
 
5611
  CharPtr pcNew;
 
5612
  Int4    i = 0, j = 0;
 
5613
  
 
5614
  pcNew = MemNew(sizeof(Char) * (Nlm_StrLen(pc) + 100));
 
5615
  
 
5616
  while (pc[i] != '\0') {
 
5617
    if (pc[i] == '>') {
 
5618
      pcNew[j] = '&';
 
5619
      pcNew[j+1] = 'g';
 
5620
      pcNew[j+2] = 't';
 
5621
      pcNew[j+3] = ';';    
 
5622
      j += 4;
 
5623
    } else if (pc[i] == '<') {
 
5624
      pcNew[j] = '&';
 
5625
      pcNew[j+1] = 'l';
 
5626
      pcNew[j+2] = 't';
 
5627
      pcNew[j+3] = ';';    
 
5628
      j += 4;
 
5629
    } else if (pc[i] == '&') {
 
5630
      pcNew[j] = '&';
 
5631
      pcNew[j+1] = 'a';
 
5632
      pcNew[j+2] = 'm';
 
5633
      pcNew[j+3] = 'p';    
 
5634
      pcNew[j+4] = ';';    
 
5635
      j += 5;
 
5636
    } else {
 
5637
      pcNew[j] = pc[i];
 
5638
      j++;    
 
5639
    }
 
5640
    i++;
 
5641
  }
 
5642
  MemFree(pc);
 
5643
  pcNew[j] = '\0';
 
5644
  pcNew = Nlm_Realloc(pcNew, sizeof(Char) * Nlm_StrLen(pcNew));
 
5645
  return(pcNew);
 
5646
}
 
5647
 
 
5648
/*---------------------------------------------------------------------------*/
 
5649
/*---------------------------------------------------------------------------*/
 
5650
/* dump out pubmed links for a CD to be used in Entrez neighboring           */
 
5651
/*---------------------------------------------------------------------------*/
 
5652
/*---------------------------------------------------------------------------*/
 
5653
void LIBCALL CddDumpPMLinks(CddPtr pcdd, FILE *FP)
 
5654
{
 
5655
  CddIdPtr     cid;
 
5656
  Int4         gi = 0, pmid = 0;
 
5657
  ValNodePtr   pub, pdescr;
 
5658
 
 
5659
  cid = (CddIdPtr) pcdd->id;
 
5660
  while (cid) {
 
5661
    switch (cid->choice) {
 
5662
      case CddId_uid:
 
5663
        gi = cid->data.intvalue;
 
5664
        break;
 
5665
      default:
 
5666
        break;
 
5667
    }
 
5668
    cid = cid->next;
 
5669
  }
 
5670
  if (gi > 0) {
 
5671
    pdescr = pcdd->description;
 
5672
    while (pdescr) {
 
5673
      switch (pdescr->choice) {
 
5674
        case CddDescr_reference:
 
5675
          pub = pdescr->data.ptrvalue;
 
5676
          if (pub->choice == PUB_PMid) {
 
5677
            pmid = pub->data.intvalue;
 
5678
            if (pmid > 0) {
 
5679
              fprintf(FP,"%d\t%d\t0\n",gi,pmid);
 
5680
            }
 
5681
          }
 
5682
          break;
 
5683
        default:
 
5684
          break;
 
5685
      }
 
5686
      pdescr = pdescr->next;
 
5687
    }
 
5688
  }
 
5689
}
 
5690
 
 
5691
/*---------------------------------------------------------------------------*/
 
5692
/*---------------------------------------------------------------------------*/
 
5693
/* dump out taxonomy links for a CD to be used in Entrez neighboring         */
 
5694
/*---------------------------------------------------------------------------*/
 
5695
/*---------------------------------------------------------------------------*/
 
5696
void LIBCALL CddDumpTaxLinks(CddPtr pcdd, FILE *FP)
 
5697
{
 
5698
  CddIdPtr     cid;
 
5699
  Int4         txid = 0, gi = 0;
 
5700
  DbtagPtr     dbtp;
 
5701
  ObjectIdPtr  oidp;
 
5702
  ValNodePtr   pdescr, vnp;
 
5703
  OrgRefPtr    pOrgRef;
 
5704
 
 
5705
  cid = (CddIdPtr) pcdd->id;
 
5706
  while (cid) {
 
5707
    switch (cid->choice) {
 
5708
      case CddId_uid:
 
5709
        gi = cid->data.intvalue;
 
5710
        break;
 
5711
      default:
 
5712
        break;
 
5713
    }
 
5714
    cid = cid->next;
 
5715
  }
 
5716
  if (gi > 0) {
 
5717
    pdescr = pcdd->description;
 
5718
    while (pdescr) {
 
5719
      switch (pdescr->choice) {
 
5720
        case CddDescr_tax_source:
 
5721
          pOrgRef = (OrgRefPtr) pdescr->data.ptrvalue;
 
5722
          vnp = pOrgRef->db;
 
5723
          while (vnp) {
 
5724
            dbtp = (DbtagPtr) vnp->data.ptrvalue;
 
5725
            if (dbtp && Nlm_StrCmp(dbtp->db,"taxon") == 0) {
 
5726
              oidp = dbtp->tag;
 
5727
              if (oidp && oidp->id >= 0) {
 
5728
                fprintf(FP,"%d\t%d\t0\n",gi,oidp->id);        
 
5729
              } 
 
5730
            }
 
5731
            vnp = vnp->next;
 
5732
          }
 
5733
          break;
 
5734
        default:
 
5735
          break;
 
5736
      }
 
5737
      pdescr = pdescr->next;
 
5738
    }
 
5739
  }
 
5740
}
 
5741
 
 
5742
/*---------------------------------------------------------------------------*/
 
5743
/*---------------------------------------------------------------------------*/
 
5744
/* dump out contents of a CD used for Entrez-indexing in pseudo-XML          */
 
5745
/*---------------------------------------------------------------------------*/
 
5746
/*---------------------------------------------------------------------------*/
 
5747
void LIBCALL CddDumpXML(CddPtr pcdd, FILE *FP)
 
5748
{
 
5749
  ValNodePtr    pdescr, cid;
 
5750
  CharPtr       cShortName  = NULL, cAbstract = NULL, cPubDate  = NULL;
 
5751
  CharPtr       cEntrezDate = NULL, cFilter   = NULL, cOrganism = NULL;
 
5752
  CharPtr       cAccession  = NULL, cDatabase = NULL;
 
5753
  Int4          gi = 0, txid = -1;
 
5754
  Boolean       bHaveAbstract = FALSE, bAllocDatabase = FALSE;
 
5755
  GlobalIdPtr   pGid;
 
5756
  OrgRefPtr     pOrgRef;
 
5757
  DatePtr       pDate;
 
5758
  ObjectIdPtr   oidp;
 
5759
  ValNodePtr    vnp;
 
5760
  DbtagPtr      dbtp;
 
5761
 
 
5762
  cEntrezDate = MemNew(11 * sizeof(Char));
 
5763
  cid = (CddIdPtr) pcdd->id;
 
5764
  while (cid) {
 
5765
    switch (cid->choice) {
 
5766
      case CddId_uid:
 
5767
        gi = cid->data.intvalue;
 
5768
        break;
 
5769
      case CddId_gid:
 
5770
        pGid = (GlobalIdPtr) cid->data.ptrvalue;
 
5771
        cAccession = pGid->accession;      
 
5772
        break;
 
5773
      default:
 
5774
        break;
 
5775
    }
 
5776
    cid = cid->next;
 
5777
  }
 
5778
  pdescr = pcdd->description;
 
5779
  while (pdescr) {
 
5780
    switch (pdescr->choice) {
 
5781
      case CddDescr_comment:
 
5782
        if (!bHaveAbstract) {
 
5783
          cAbstract = (CharPtr) pdescr->data.ptrvalue;
 
5784
          bHaveAbstract = TRUE;
 
5785
        }
 
5786
        break;
 
5787
      case CddDescr_source:
 
5788
        cDatabase = (CharPtr) pdescr->data.ptrvalue;
 
5789
      case CddDescr_tax_source:
 
5790
        pOrgRef = (OrgRefPtr) pdescr->data.ptrvalue;
 
5791
        cOrganism = pOrgRef->taxname; 
 
5792
        vnp = pOrgRef->db;
 
5793
        while (vnp) {
 
5794
          dbtp = (DbtagPtr) vnp->data.ptrvalue;
 
5795
          if (dbtp && Nlm_StrCmp(dbtp->db,"taxon") == 0) {
 
5796
            oidp = dbtp->tag;
 
5797
            if (oidp && oidp->id >= 0) {
 
5798
              txid = oidp->id;        
 
5799
            } 
 
5800
          }
 
5801
          vnp = vnp->next;
 
5802
        }
 
5803
        break;
 
5804
      case CddDescr_create_date:
 
5805
      case CddDescr_update_date:
 
5806
        cPubDate = MemNew(11 * sizeof(Char));
 
5807
        pDate = (DatePtr) pdescr->data.ptrvalue;
 
5808
        if (pDate->data[2] > 0 && pDate->data[3] > 0)
 
5809
          sprintf(cPubDate,"%d/%d/%d",pDate->data[1]+1900,pDate->data[2],pDate->data[3]);
 
5810
        if (pDate->data[2] > 0 && pDate->data[3] == 0)
 
5811
          sprintf(cPubDate,"%d/%d",pDate->data[1]+1900,pDate->data[2]);
 
5812
        if (pDate->data[2] == 0 && pDate->data[3] == 0)
 
5813
          sprintf(cPubDate,"%d",pDate->data[1]+1900);
 
5814
        break;
 
5815
      default:
 
5816
        break;
 
5817
    }
 
5818
    pdescr = pdescr->next;
 
5819
  }
 
5820
  cShortName = (CharPtr) pcdd->name;
 
5821
  if (cAccession) {
 
5822
    if (Nlm_StrNCmp(cAccession,"cd",2) == 0) {
 
5823
      cDatabase = StringSave("Cdd");
 
5824
      bAllocDatabase = TRUE;
 
5825
    }
 
5826
  }
 
5827
  if (!cPubDate) {
 
5828
    cPubDate = MemNew(11 * sizeof(Char));
 
5829
    pDate = (DatePtr) DateCurr();
 
5830
    if (pDate->data[2] > 0 && pDate->data[3] > 0)
 
5831
      sprintf(cPubDate,"%d/%d/%d",pDate->data[1]+1900,pDate->data[2],pDate->data[3]);
 
5832
    if (pDate->data[2] > 0 && pDate->data[3] == 0)
 
5833
      sprintf(cPubDate,"%d/%d",pDate->data[1]+1900,pDate->data[2]);
 
5834
    if (pDate->data[2] == 0 && pDate->data[3] == 0)
 
5835
      sprintf(cPubDate,"%d",pDate->data[1]+1900);
 
5836
  }
 
5837
  pDate = (DatePtr) DateCurr();
 
5838
  if (pDate->data[2] > 0 && pDate->data[3] > 0)
 
5839
    sprintf(cEntrezDate,"%d/%d/%d",pDate->data[1]+1900,pDate->data[2],pDate->data[3]);
 
5840
  if (pDate->data[2] > 0 && pDate->data[3] == 0)
 
5841
    sprintf(cEntrezDate,"%d/%d",pDate->data[1]+1900,pDate->data[2]);
 
5842
  if (pDate->data[2] == 0 && pDate->data[3] == 0)
 
5843
    sprintf(cEntrezDate,"%d",pDate->data[1]+1900);
 
5844
 
 
5845
  cAbstract = CddFixForXML(cAbstract);
 
5846
  cShortName = CddFixForXML(cShortName);
 
5847
 
 
5848
  fprintf(FP,"    <IdxDocument>\n");  
 
5849
  if (gi > 0) 
 
5850
    fprintf(FP,"        <IdxUid>%d</IdxUid>\n",gi);
 
5851
  fprintf(FP,"        <IdxKeyFields>\n");
 
5852
  if (cPubDate)
 
5853
    fprintf(FP,"            <PubDate>%s</PubDate>\n",cPubDate);  
 
5854
  if (cEntrezDate)
 
5855
    fprintf(FP,"            <EntrezDate>%s</EntrezDate>\n",cEntrezDate);  
 
5856
  fprintf(FP,"        </IdxKeyFields>\n");
 
5857
  fprintf(FP,"        <IdxDisplayFields>\n");
 
5858
  if (cAccession)
 
5859
    fprintf(FP,"            <Accession>%s</Accession>\n",cAccession);
 
5860
  if (cShortName)
 
5861
    fprintf(FP,"            <Title>%s</Title>\n",cShortName);
 
5862
  if (cAbstract) 
 
5863
    fprintf(FP,"            <Abstract>%s</Abstract>\n",cAbstract);
 
5864
  fprintf(FP,"        </IdxDisplayFields>\n");
 
5865
  fprintf(FP,"        <IdxSearchFields>\n");
 
5866
  if (gi > 0) 
 
5867
    fprintf(FP,"            <Uid>%d</Uid>\n",gi);
 
5868
  if (cFilter)
 
5869
    fprintf(FP,"            <Filter>%s</Filter>\n",cFilter);
 
5870
  if (cAccession)
 
5871
    fprintf(FP,"            <Accession>%s</Accession>\n",cAccession);
 
5872
  if (cDatabase) 
 
5873
    fprintf(FP,"            <Database>%s</Database>\n",cDatabase);
 
5874
  if (cShortName)
 
5875
    fprintf(FP,"            <Title>%s</Title>\n",cShortName);
 
5876
  if (cAbstract) 
 
5877
    fprintf(FP,"            <Text>%s</Text>\n",cAbstract);
 
5878
  if (txid >= 0) 
 
5879
    fprintf(FP,"            <Organism>txid%d</Organism>\n",txid);
 
5880
  if (cOrganism)
 
5881
    fprintf(FP,"            <Organism>%s</Organism>\n",cOrganism);
 
5882
  if (cPubDate)
 
5883
    fprintf(FP,"            <PubDate>%s</PubDate>\n",cPubDate);  
 
5884
  if (cEntrezDate)
 
5885
    fprintf(FP,"            <EntrezDate>%s</EntrezDate>\n",cEntrezDate);  
 
5886
  fprintf(FP,"        </IdxSearchFields>\n");
 
5887
  fprintf(FP,"    </IdxDocument>\n");  
 
5888
 
 
5889
  if (bAllocDatabase) MemFree(cDatabase);
 
5890
  MemFree(cPubDate);
 
5891
  MemFree(cEntrezDate);
 
5892
}
 
5893
 
 
5894
/*---------------------------------------------------------------------------*/
 
5895
/*---------------------------------------------------------------------------*/
 
5896
/* create and link a CddIdxData structure                                    */
 
5897
/*---------------------------------------------------------------------------*/
 
5898
/*---------------------------------------------------------------------------*/
 
5899
CddIdxDataPtr LIBCALL CddIdxDataNew(CharPtr acc, Int4 uid)
 
5900
{
 
5901
 
 
5902
  CddIdxDataPtr     cidp;
 
5903
 
 
5904
  cidp = MemNew(sizeof(CddIdxData));
 
5905
 
 
5906
  if (acc) {
 
5907
    cidp->cCDDid  = acc;
 
5908
    cidp->iPssmId = uid;
 
5909
  }
 
5910
  cidp->next = NULL;
 
5911
  return(cidp);
 
5912
}
 
5913
 
 
5914
CddIdxDataPtr LIBCALL CddIdxDataLink(CddIdxDataPtr PNTR head, CddIdxDataPtr cidp)
 
5915
{
 
5916
  CddIdxDataPtr  cidpThis;
 
5917
 
 
5918
  cidpThis = *head;
 
5919
  if (cidpThis) {
 
5920
    while (cidpThis->next != NULL) cidpThis = cidpThis->next;
 
5921
    cidpThis->next = cidp;
 
5922
  } else {
 
5923
    *head = cidp;
 
5924
  }
 
5925
  return *head;
 
5926
}
 
5927
 
 
5928
 
 
5929
/*---------------------------------------------------------------------------*/
 
5930
/*---------------------------------------------------------------------------*/
 
5931
/* read in the cdd.idx file to return data associating accessions with uids  */
 
5932
/*---------------------------------------------------------------------------*/
 
5933
/*---------------------------------------------------------------------------*/
 
5934
CddIdxDataPtr LIBCALL CddReadIdx(CharPtr CDDidx)
 
5935
{
 
5936
  CddIdxDataPtr cidp, cidpHead = NULL;
 
5937
  Char          pcBuf[2048];
 
5938
  CharPtr       pcTest;
 
5939
  CharPtr       acc;
 
5940
  Int4          uid;
 
5941
  FILE         *fp;
 
5942
  
 
5943
  fp = FileOpen(CDDidx,"r");
 
5944
  if (!fp) CddSevError("could not open cdd.idx!");
 
5945
  do {
 
5946
    pcBuf[0]='\0';
 
5947
    pcTest = fgets(pcBuf, (size_t)2048, fp);
 
5948
    if (pcTest) if (pcTest[0] != ' ') {
 
5949
      uid = (Int4) atoi(Nlm_StrTok(pcTest," "));
 
5950
      acc = (CharPtr) StringSave(Nlm_StrTok(NULL," "));
 
5951
      if (acc) {
 
5952
        acc[Nlm_StrLen(acc) - 1] = '\0';
 
5953
        cidp = CddIdxDataNew(acc,uid);
 
5954
        CddIdxDataLink(&(cidpHead),cidp);
 
5955
      }
 
5956
    }
 
5957
  } while (pcTest);
 
5958
  FileClose(fp);
 
5959
  return(cidpHead);
 
5960
}
 
5961
 
 
5962
 
 
5963
/*---------------------------------------------------------------------------*/
 
5964
/*---------------------------------------------------------------------------*/
 
5965
/* retrieve a CDD accession via the PSSMid or vice versa                     */
 
5966
/*---------------------------------------------------------------------------*/
 
5967
/*---------------------------------------------------------------------------*/
 
5968
void LIBCALL CddAccFromPssmId(Int4 iPssmId, CharPtr cAcc, CharPtr CDDidx)
 
5969
{
 
5970
  CddIdxDataPtr cidp;
 
5971
 
 
5972
  cidp = CddReadIdx(CDDidx);
 
5973
  if (!cidp) CddSevError("cdd.idx read failed!");
 
5974
  while (cidp) {
 
5975
    if (iPssmId == cidp->iPssmId) {
 
5976
      Nlm_StrCpy(cAcc, cidp->cCDDid);
 
5977
      break;
 
5978
    }
 
5979
    cidp = cidp->next;
 
5980
  }
 
5981
}
 
5982
 
 
5983
void LIBCALL CddPssmIdFromAcc(Int4 *iPssmId, CharPtr cAcc, CharPtr CDDidx)
 
5984
{
 
5985
  CddIdxDataPtr cidp;
 
5986
 
 
5987
  cidp = CddReadIdx(CDDidx);
 
5988
  if (!cidp) CddSevError("cdd.idx read failed!");
 
5989
  while (cidp) {
 
5990
    if (Nlm_StrCmp(cAcc, cidp->cCDDid) == 0) {
 
5991
      *iPssmId = cidp->iPssmId;
 
5992
      break;
 
5993
    }
 
5994
    cidp = cidp->next;
 
5995
  }
 
5996
}
 
5997
 
 
5998
/*---------------------------------------------------------------------------*/
 
5999
/*---------------------------------------------------------------------------*/
 
6000
/* truncate a string right at the position of the first puncutation mark     */
 
6001
/*---------------------------------------------------------------------------*/
 
6002
/*---------------------------------------------------------------------------*/
 
6003
void LIBCALL CddTruncStringAtFirstPunct(CharPtr pChar)
 
6004
{
 
6005
  Int4   i;
 
6006
  
 
6007
  for (i=0;i<Nlm_StrLen(pChar);i++) {
 
6008
    if (pChar[i] == ',' ||
 
6009
        pChar[i] == ';' ||
 
6010
        pChar[i] == '.') {
 
6011
      pChar[i] = '\0'; break;   
 
6012
    }
 
6013
  }
 
6014
}
 
6015
 
 
6016
void LIBCALL CddFillBlanksInString(CharPtr pChar) 
 
6017
{
 
6018
  Int4   i;
 
6019
  
 
6020
  for (i=0;i<Nlm_StrLen(pChar);i++) {
 
6021
    if (pChar[i] == ' ') pChar[i] = '_';
 
6022
  }
 
6023
}
 
6024
 
 
6025
/*---------------------------------------------------------------------------*/
 
6026
/*---------------------------------------------------------------------------*/
 
6027
/* Remove a Consensus from a CD                                              */
 
6028
/*---------------------------------------------------------------------------*/
 
6029
/*---------------------------------------------------------------------------*/
 
6030
Boolean LIBCALL CddRemoveConsensus(CddPtr pcdd)
 
6031
{
 
6032
  SeqIdPtr     sip;
 
6033
  BioseqPtr    bsp;
 
6034
  SeqAlignPtr  salp, salpnew;
 
6035
  DenseDiagPtr ddp;
 
6036
  BioseqSetPtr bssp;
 
6037
  SeqEntryPtr  sepThis, sepHead = NULL, sepTail = NULL;
 
6038
 
 
6039
  if (!CddHasConsensus(pcdd)) return FALSE;
 
6040
  salp = (SeqAlignPtr)  pcdd->seqannot->data;
 
6041
  ddp  = (DenseDiagPtr) salp->segs;
 
6042
  sip  = (SeqIdPtr)     ddp->id->next;
 
6043
  bssp = (BioseqSetPtr) pcdd->sequences->data.ptrvalue;
 
6044
  if (pcdd->alignannot) {
 
6045
    CddTransferAlignAnnot(pcdd->alignannot, sip, salp, bssp);
 
6046
  }
 
6047
  salpnew = CddReindexSeqAlign(salp,sip,bssp);
 
6048
  pcdd->seqannot->data = salpnew->next;
 
6049
  sepThis = bssp->seq_set;
 
6050
  while (sepThis) {
 
6051
    bsp = (BioseqPtr) sepThis->data.ptrvalue;
 
6052
    if (!SipIsConsensus(bsp->id)) {
 
6053
      if (!sepHead) sepHead = sepThis;
 
6054
      if (sepTail) sepTail->next = sepThis;
 
6055
      sepTail = sepThis;
 
6056
    }
 
6057
    sepThis = sepThis->next;
 
6058
  }
 
6059
  sepTail->next = NULL;
 
6060
  bssp->seq_set = sepHead;
 
6061
  pcdd->sequences->data.ptrvalue = bssp;
 
6062
  pcdd->profile_range = NULL;
 
6063
  pcdd->trunc_master = NULL;
 
6064
  pcdd->posfreq = NULL;
 
6065
  pcdd->scoremat = NULL;
 
6066
  pcdd->distance = NULL;
 
6067
  return TRUE;
 
6068
}
 
6069
 
 
6070
/*---------------------------------------------------------------------------*/
 
6071
/*---------------------------------------------------------------------------*/
 
6072
/* Check if a CD has annotation that is inconsistent with the alignment      */
 
6073
/*---------------------------------------------------------------------------*/
 
6074
/*---------------------------------------------------------------------------*/
 
6075
Boolean LIBCALL CddFeaturesAreConsistent(CddPtr pcdd, CharPtr errmsg)
 
6076
{
 
6077
  SeqAlignPtr    salp;
 
6078
  CddExpAlignPtr eap;
 
6079
  BioseqSetPtr   bssp;
 
6080
  AlignAnnotPtr  aap;
 
6081
  SeqLocPtr      slp;
 
6082
  Int4           iWhere;
 
6083
  
 
6084
  if (!pcdd->alignannot) return TRUE;
 
6085
  bssp = (BioseqSetPtr) pcdd->sequences->data.ptrvalue;
 
6086
  salp = (SeqAlignPtr) pcdd->seqannot->data;
 
6087
  if (!salp) return TRUE;
 
6088
  eap = SeqAlignToCddExpAlign(salp,bssp->seq_set);
 
6089
  aap = pcdd->alignannot;
 
6090
  while (aap) {
 
6091
    slp = aap->location;
 
6092
    iWhere = CddSeqLocInExpAlign(slp,eap);
 
6093
    if (iWhere >= 0) {
 
6094
      if (iWhere == INT4_MAX) {
 
6095
        sprintf(errmsg,"Feature: %s, illegal feature address",aap->description);
 
6096
      } else {
 
6097
        sprintf(errmsg,"Feature: %s, Location: %d",aap->description,iWhere+1);
 
6098
      }    
 
6099
      return FALSE;
 
6100
    }
 
6101
    aap = aap->next;
 
6102
  }
 
6103
  CddExpAlignFree(eap);
 
6104
  return TRUE;
 
6105
  
 
6106
}
 
6107
 
 
6108
/*---------------------------------------------------------------------------*/
 
6109
/*---------------------------------------------------------------------------*/
 
6110
/* find a MMDB-Identifier in a Bioseq                                        */
 
6111
/*---------------------------------------------------------------------------*/
 
6112
/*---------------------------------------------------------------------------*/
 
6113
SeqAnnotPtr LIBCALL CddFindMMDBIdInBioseq(BioseqPtr bsp, Int4 *iMMDBid)
 
6114
{
 
6115
  SeqAnnotPtr annot, prevannot = NULL;
 
6116
  SeqIdPtr    sip;
 
6117
  DbtagPtr    dbtp;
 
6118
  ObjectIdPtr oidp;
 
6119
 
 
6120
  *iMMDBid = 0;
 
6121
  annot = bsp->annot;
 
6122
  while(annot) {
 
6123
    if (annot->type == 4) {
 
6124
      sip = annot->data;
 
6125
      while (sip) {
 
6126
        if (sip->choice == SEQID_GENERAL) {
 
6127
          dbtp = sip->data.ptrvalue;
 
6128
          if (Nlm_StrCmp(dbtp->db,"mmdb") == 0) {
 
6129
            oidp = dbtp->tag;
 
6130
            *iMMDBid = oidp->id;
 
6131
            return(prevannot);
 
6132
          }
 
6133
        }
 
6134
        sip = sip->next;
 
6135
      }
 
6136
    }
 
6137
    if (annot->next) prevannot = annot;
 
6138
    annot = annot->next;
 
6139
  }
 
6140
  return (NULL);
 
6141
}
 
6142
 
 
6143
 
 
6144
/*---------------------------------------------------------------------------*/
 
6145
/*---------------------------------------------------------------------------*/
 
6146
/* Check if a CD has 3D superposition information for all aligned chains     */
 
6147
/*---------------------------------------------------------------------------*/
 
6148
/*---------------------------------------------------------------------------*/
 
6149
Boolean LIBCALL CddHas3DSuperpos(CddPtr pcdd)
 
6150
{
 
6151
  BioseqSetPtr          bssp;
 
6152
  BiostrucAnnotSetPtr   basp        = NULL;
 
6153
  SeqAlignPtr           salp;
 
6154
  DenseDiagPtr          ddp;
 
6155
  SeqIdPtr              sip1, sip2;
 
6156
  Int4Ptr               pMMDBid;
 
6157
  Int4                  i, nStruct  = 0;
 
6158
  BioseqPtr             bsp;
 
6159
  BiostrucFeatureSetPtr bfsp;
 
6160
  BiostrucFeaturePtr    bfp;
 
6161
  ValNodePtr            location, vnp;
 
6162
  ChemGraphAlignmentPtr cgap;
 
6163
  Int4                  thisslave;
 
6164
 
 
6165
  pMMDBid = MemNew(250 * sizeof(Int4));
 
6166
  basp = pcdd->features;
 
6167
  if (!basp) return FALSE;
 
6168
  bssp = (BioseqSetPtr) pcdd->sequences->data.ptrvalue;
 
6169
/*---------------------------------------------------------------------------*/
 
6170
/* create a list of MMDB-Id pairs which need to be checked                   */
 
6171
/*---------------------------------------------------------------------------*/
 
6172
  salp = pcdd->seqannot->data;
 
6173
  while (salp) {
 
6174
    ddp = salp->segs;
 
6175
    sip1 = ddp->id;
 
6176
          sip2 = ddp->id->next;
 
6177
    if (sip1->choice == SEQID_PDB) {
 
6178
                  if (nStruct == 0) {
 
6179
        bsp = CddFindSip(sip1,bssp->seq_set);
 
6180
        CddFindMMDBIdInBioseq(bsp,&pMMDBid[nStruct]);
 
6181
        nStruct++;
 
6182
      }
 
6183
    }
 
6184
    if (sip2->choice == SEQID_PDB) {
 
6185
      bsp = CddFindSip(sip2,bssp->seq_set);
 
6186
      CddFindMMDBIdInBioseq(bsp,&pMMDBid[nStruct]);
 
6187
      nStruct++;
 
6188
    }
 
6189
    salp = salp->next;
 
6190
  }
 
6191
/*---------------------------------------------------------------------------*/
 
6192
/* now walk down the biostruc annot set and remove all pairs encountered     */
 
6193
/*---------------------------------------------------------------------------*/
 
6194
  bfsp = basp->features;
 
6195
  while (bfsp) {
 
6196
    bfp = bfsp->features;
 
6197
    while (bfp) {
 
6198
      location = bfp->Location_location;
 
6199
      if (location->choice == Location_location_alignment) {
 
6200
        cgap = location->data.ptrvalue;
 
6201
        vnp = cgap->biostruc_ids;
 
6202
        thisslave = vnp->next->data.intvalue;
 
6203
        for (i=1;i<nStruct;i++) {
 
6204
          if (thisslave == pMMDBid[i]) {
 
6205
            pMMDBid[i] = 0; break;  
 
6206
          }
 
6207
        }
 
6208
      }
 
6209
      bfp = bfp->next;    
 
6210
    }
 
6211
    bfsp = bfsp->next;  
 
6212
  }
 
6213
 
 
6214
/*---------------------------------------------------------------------------*/
 
6215
/* at this point check if any pair is left (i.e. any mmdb-id except master)  */
 
6216
/*---------------------------------------------------------------------------*/
 
6217
  for (i=1;i<nStruct;i++) {
 
6218
    if (pMMDBid[i]) {
 
6219
      MemFree(pMMDBid); return(FALSE);
 
6220
    }
 
6221
  }
 
6222
  MemFree(pMMDBid); return TRUE;
 
6223
}
 
6224
 
 
6225
/*---------------------------------------------------------------------------*/
 
6226
/*---------------------------------------------------------------------------*/
 
6227
/* CD has unfinished alignment business                                      */
 
6228
/*---------------------------------------------------------------------------*/
 
6229
/*---------------------------------------------------------------------------*/
 
6230
Boolean    LIBCALL CddHasPendingAlignments(CddPtr pcdd)
 
6231
{
 
6232
  if (pcdd->pending) return TRUE;
 
6233
  return FALSE;
 
6234
}
 
6235
 
 
6236
/*---------------------------------------------------------------------------*/
 
6237
/*---------------------------------------------------------------------------*/
 
6238
/* check 
 
6239
/*---------------------------------------------------------------------------*/
 
6240
/*---------------------------------------------------------------------------*/
 
6241
Boolean LIBCALL SeqHasTax(BioseqPtr bsp)
 
6242
{
 
6243
  ValNodePtr  descr;
 
6244
 
 
6245
  descr = bsp->descr;
 
6246
  while (descr) {
 
6247
    if (descr->choice == Seq_descr_source) {
 
6248
      return(TRUE);
 
6249
    }
 
6250
    descr = descr->next;
 
6251
  }
 
6252
  return(FALSE);
 
6253
}
 
6254
 
 
6255
/*---------------------------------------------------------------------------*/
 
6256
/*---------------------------------------------------------------------------*/
 
6257
/* Fetch bioseq using BLAST dB, removing ids redundant to query.             */
 
6258
/* index is the NR index in the dB -- if you don't have it, use -1.          */
 
6259
/*---------------------------------------------------------------------------*/
 
6260
/*---------------------------------------------------------------------------*/
 
6261
NLM_EXTERN BlastDefLineSetPtr LIBCALL BlastDefLineSetFree PROTO ((BlastDefLineSetPtr ));
 
6262
 
 
6263
BioseqPtr LIBCALL CddReadDBGetBioseq(SeqIdPtr query, Int4 index,
 
6264
                                     ReadDBFILEPtr rdfp)
 
6265
{
 
6266
  return CddReadDBGetBioseqEx(query,index,rdfp,TRUE);
 
6267
}
 
6268
 
 
6269
BioseqPtr LIBCALL CddReadDBGetBioseqEx(SeqIdPtr query, Int4 index,
 
6270
                                       ReadDBFILEPtr rdfp, Boolean bUseObjMgr)
 
6271
{
 
6272
    BioseqPtr           bsp;
 
6273
    BlastDefLinePtr     bdp,tbdp;
 
6274
    RDBTaxNamesPtr      tnames;
 
6275
    Char                buf[64];
 
6276
 
 
6277
    if(query == NULL) return NULL;
 
6278
    if(index < 0 && (index = SeqId2OrdinalId(rdfp, query)) < 0)
 
6279
        return NULL;
 
6280
    if(!(bsp = readdb_get_bioseq_ex(rdfp, index, bUseObjMgr, FALSE))) return NULL;
 
6281
    if(!(tbdp = FDGetDeflineAsnFromBioseq(bsp))) return NULL;
 
6282
 
 
6283
    for(bdp=tbdp; bdp; bdp=bdp->next) {
 
6284
        if(SeqIdComp(bdp->seqid,query) == SIC_YES ||
 
6285
           SeqIdComp(bdp->seqid->next,query) == SIC_YES) {
 
6286
            tnames = RDBGetTaxNames(rdfp->taxinfo, bdp->taxid);
 
6287
            break;
 
6288
        }
 
6289
    }
 
6290
    if(bdp == NULL) {
 
6291
        SeqIdPrint(query,buf,PRINTID_FASTA_SHORT);
 
6292
        ErrPostEx(SEV_WARNING,0,0,"Bad seqid %s",buf);
 
6293
        return NULL;
 
6294
    }
 
6295
    bsp->id = SeqIdSetFree(bsp->id);
 
6296
    bsp->id = SeqIdDup(bdp->seqid->next);
 
6297
    bsp->id->next = SeqIdDup(bdp->seqid);
 
6298
 
 
6299
    /* Remove concatenated deflines. */
 
6300
    SeqDescrFree(ValNodeExtractList(&bsp->descr, Seq_descr_title));
 
6301
    /* Replace with single defline. */
 
6302
    ValNodeCopyStr(&bsp->descr, Seq_descr_title, bdp->title);
 
6303
 
 
6304
    tbdp = BlastDefLineSetFree(tbdp);
 
6305
 
 
6306
    /* Remove all user data (tax info, asn1 defline). */
 
6307
    SeqDescrFree(ValNodeExtractList(&bsp->descr, Seq_descr_user));
 
6308
 
 
6309
    if(tnames) {                /* Add taxid, taxnames. */
 
6310
        BioSourcePtr    bsrp;
 
6311
        DbtagPtr        dtp;
 
6312
 
 
6313
        bsrp = BioSourceNew();
 
6314
        bsrp->subtype = NULL;
 
6315
        bsrp->org = OrgRefNew();
 
6316
        if(tnames->sci_name)
 
6317
            bsrp->org->taxname = StringSave(tnames->sci_name);
 
6318
        if(tnames->common_name)
 
6319
            bsrp->org->common = StringSave(tnames->common_name);
 
6320
        dtp = DbtagNew();
 
6321
        dtp->db = StringSave("taxon");
 
6322
        dtp->tag = ObjectIdNew();
 
6323
        dtp->tag->id = tnames->tax_id;
 
6324
        bsrp->org->db = NULL;
 
6325
        ValNodeAddPointer(&bsrp->org->db, 0, dtp);
 
6326
        ValNodeAddPointer(&bsp->descr, Seq_descr_source, bsrp);
 
6327
    }
 
6328
    RDBTaxNamesFree(tnames);
 
6329
 
 
6330
/*    if (!SeqHasTax(bsp)) {
 
6331
      printf("Warning: Bioseq without taxonomy!");
 
6332
      if (query->choice == SEQID_GI) {
 
6333
        printf("GI: %d\n",query->data.intvalue);
 
6334
      } else printf("\n");
 
6335
    } */
 
6336
    return bsp;
 
6337
}
 
6338
 
 
6339
/*---------------------------------------------------------------------------*/
 
6340
/*---------------------------------------------------------------------------*/
 
6341
/* return a common style dictionary for Cn3D 4.0                             */
 
6342
/*---------------------------------------------------------------------------*/
 
6343
/*---------------------------------------------------------------------------*/
 
6344
/*
 
6345
Cn3dStyleDictionaryPtr LIBCALL CddSrvGetStyle2(Int4 *styles[], Int4 nstyles)
 
6346
{
 
6347
    Cn3dStyleDictionaryPtr    csdp;
 
6348
    Cn3dStyleTableItemPtr     cstip,cstipTail,cstipHead=NULL;
 
6349
    Int4                            i;
 
6350
 
 
6351
    if(nstyles < 1) return(NULL);
 
6352
    csdp = Cn3dStyleDictionaryNew();
 
6353
    csdp->global_style = CddSrvGetStyle2_Ex(styles[0]);
 
6354
    for(i=1; i<nstyles; i++) {
 
6355
        cstip = Cn3dStyleTableItemNew();
 
6356
        cstip->id = i;
 
6357
        cstip->style = CddSrvGetStyle2_Ex(styles[i]);
 
6358
        cstip->next = NULL;
 
6359
        if(!cstipHead) cstipHead = cstipTail = cstip;
 
6360
        else {
 
6361
            cstipTail->next = cstip;
 
6362
            cstipTail = cstip;
 
6363
        }
 
6364
    }
 
6365
    csdp->style_table = cstipHead;
 
6366
 
 
6367
    return(csdp);
 
6368
}
 
6369
static Cn3dStyleSettingsPtr CddSrvGetStyle2_Ex(Int4 style[])
 
6370
{
 
6371
  Cn3dStyleSettingsPtr      cssp;
 
6372
  Cn3dBackboneStylePtr      cbsp;
 
6373
  Cn3dGeneralStylePtr       cgsp;
 
6374
  Cn3dBackboneLabelStylePtr cblsp;
 
6375
 
 
6376
    cssp = Cn3dStyleSettingsNew();
 
6377
      cssp->next = NULL;
 
6378
      cssp->name = NULL;
 
6379
      cbsp = Cn3dBackboneStyleNew();
 
6380
        cbsp->type = style[0];
 
6381
        cbsp->style = style[1];
 
6382
        cbsp->color_scheme = style[2];
 
6383
        cbsp->user_color =
 
6384
            MyCn3dColorInit(style[3],style[4],style[5],
 
6385
                            style[6],style[7]);
 
6386
      cssp->protein_backbone = cbsp;
 
6387
      cbsp = Cn3dBackboneStyleNew();
 
6388
        cbsp->type = style[8];
 
6389
        cbsp->style = style[9];
 
6390
        cbsp->color_scheme = style[10];
 
6391
        cbsp->user_color =
 
6392
            MyCn3dColorInit(style[11],style[12],style[13],
 
6393
                            style[14],style[15]);
 
6394
      cssp->nucleotide_backbone = cbsp;
 
6395
        cgsp = Cn3dGeneralStyleNew();
 
6396
        cgsp->is_on = style[16];
 
6397
        cgsp->style = style[17];
 
6398
        cgsp->color_scheme = style[18];
 
6399
        cgsp->user_color =
 
6400
            MyCn3dColorInit(style[19],style[20],style[21],
 
6401
                            style[22],style[23]);
 
6402
      cssp->protein_sidechains = cgsp;
 
6403
        cgsp = Cn3dGeneralStyleNew();
 
6404
        cgsp->is_on = style[24];
 
6405
        cgsp->style = style[25];
 
6406
        cgsp->color_scheme = style[26];
 
6407
        cgsp->user_color =
 
6408
            MyCn3dColorInit(style[27],style[28],style[29],
 
6409
                            style[30],style[31]);
 
6410
      cssp->nucleotide_sidechains = cgsp;
 
6411
        cgsp = Cn3dGeneralStyleNew();
 
6412
        cgsp->is_on = style[32];
 
6413
        cgsp->style = style[33];
 
6414
        cgsp->color_scheme = style[34];
 
6415
        cgsp->user_color =
 
6416
            MyCn3dColorInit(style[35],style[36],style[37],
 
6417
                            style[38],style[39]);
 
6418
       cssp->heterogens = cgsp;
 
6419
        cgsp = Cn3dGeneralStyleNew();
 
6420
        cgsp->is_on = style[40];
 
6421
        cgsp->style = style[41];
 
6422
        cgsp->color_scheme = style[42];
 
6423
        cgsp->user_color =
 
6424
            MyCn3dColorInit(style[43],style[44],style[45],
 
6425
                            style[46],style[47]);
 
6426
      cssp->solvents = cgsp;
 
6427
        cgsp = Cn3dGeneralStyleNew();
 
6428
        cgsp->is_on = style[48];
 
6429
        cgsp->style = style[49];
 
6430
        cgsp->color_scheme = style[50];
 
6431
        cgsp->user_color =
 
6432
            MyCn3dColorInit(style[51],style[52],style[53],
 
6433
                            style[54],style[55]);
 
6434
       cssp->connections = cgsp;
 
6435
        cgsp = Cn3dGeneralStyleNew();
 
6436
        cgsp->is_on = style[56];
 
6437
        cgsp->style = style[57];
 
6438
        cgsp->color_scheme = style[58];
 
6439
        cgsp->user_color =
 
6440
            MyCn3dColorInit(style[59],style[60],style[61],
 
6441
                            style[62],style[63]);
 
6442
      cssp->helix_objects = cgsp;
 
6443
        cgsp = Cn3dGeneralStyleNew();
 
6444
        cgsp->is_on = style[64];
 
6445
        cgsp->style = style[65];
 
6446
        cgsp->color_scheme = style[66];
 
6447
        cgsp->user_color =
 
6448
            MyCn3dColorInit(style[67],style[68],style[69],
 
6449
                            style[70],style[71]);
 
6450
      cssp->strand_objects = cgsp;
 
6451
      cssp->virtual_disulfides_on = style[72];
 
6452
      cssp->virtual_disulfide_color =
 
6453
          MyCn3dColorInit(style[73],style[74],style[75],
 
6454
                          style[76],style[77]);
 
6455
      cssp->hydrogens_on = style[78];
 
6456
      cssp->background_color =
 
6457
          MyCn3dColorInit(style[79],style[80],style[81],
 
6458
                          style[82],style[83]);
 
6459
      cssp->scale_factor = style[84];
 
6460
      cssp->space_fill_proportion = style[85];
 
6461
      cssp->ball_radius = style[86];
 
6462
      cssp->stick_radius = style[87];
 
6463
      cssp->tube_radius = style[88];
 
6464
      cssp->tube_worm_radius = style[89];
 
6465
      cssp->helix_radius = style[90];
 
6466
      cssp->strand_width = style[91];
 
6467
      cssp->strand_thickness = style[92];
 
6468
        cblsp = Cn3dBackboneLabelStyleNew();
 
6469
        cblsp->spacing = style[93];
 
6470
        cblsp->type = style[94];
 
6471
        cblsp->number = style[95];
 
6472
        cblsp->termini = style[96];
 
6473
        cblsp->white = style[97];
 
6474
      cssp->protein_labels = cblsp;
 
6475
        cblsp = Cn3dBackboneLabelStyleNew();
 
6476
        cblsp->spacing = style[98];
 
6477
        cblsp->type = style[99];
 
6478
        cblsp->number = style[100];
 
6479
        cblsp->termini = style[101];
 
6480
        cblsp->white = style[102];
 
6481
      cssp->nucleotide_labels =cblsp;
 
6482
      cssp->ion_labels = style[103];
 
6483
 
 
6484
  return(cssp);
 
6485
}
 
6486
static Cn3dColorPtr MyCn3dColorInit(Int4 scale_factor, Int4 red, Int4 green, Int4 blue,
 
6487
                             Int4 alpha)
 
6488
{
 
6489
    Cn3dColorPtr        ccp;
 
6490
 
 
6491
    ccp = Cn3dColorNew();
 
6492
    ccp->scale_factor = scale_factor;
 
6493
    ccp->red = red;
 
6494
    ccp->green = green;
 
6495
    ccp->blue = blue;
 
6496
    ccp->alpha = alpha;
 
6497
 
 
6498
    return ccp;
 
6499
}
 
6500
*/