1 static char const rcsid[] = "$Id: urkbias.c,v 6.7 2003/05/30 17:25:38 coulouri Exp $";
2
3 /*
4 * ===========================================================================
5 *
6 * PUBLIC DOMAIN NOTICE
7 * National Center for Biotechnology Information
8 *
9 * This software/database is a "United States Government Work" under the
10 * terms of the United States Copyright Act. It was written as part of
11 * the author's official duties as a United States Government employee and
12 * thus cannot be copyrighted. This software/database is freely available
13 * to the public for use. The National Library of Medicine and the U.S.
14 * Government have not placed any restriction on its use or reproduction.
15 *
16 * Although all reasonable efforts have been taken to ensure the accuracy
17 * and reliability of the software and data, the NLM and the U.S.
18 * Government do not and cannot warrant the performance or results that
19 * may be obtained by using this software or data. The NLM and the U.S.
20 * Government disclaim all warranties, express or implied, including
21 * warranties of performance, merchantability or fitness for any particular
22 * purpose.
23 *
24 * Please cite the author in any work or product based on this material.
25 *
26 * ===========================================================================
27 *
28 * File Name: urkbias.c
29 *
30 * Author(s): John Kuzio
31 *
32 * Version Creation Date: 98-01-01
33 *
34 * $Revision: 6.7 $
35 *
36 * File Description: codon bias
37 *
38 * Modifications:
39 * --------------------------------------------------------------------------
40 * Date Name Description of modification
41 * --------------------------------------------------------------------------
42 * $Log: urkbias.c,v $
43 * Revision 6.7 2003/05/30 17:25:38 coulouri
44 * add rcsid
45 *
46 * Revision 6.6 1998/09/30 21:46:45 kuzio
47 * no need to reverse bias score
48 *
49 * Revision 6.5 1998/09/16 18:03:31 kuzio
50 * cvs logging
51 *
52 *
53 * ==========================================================================
54 */
55
56 #include <ncbi.h>
57 #include <accentr.h>
58 #include <gather.h>
59 #include <urkcnsrt.h>
60 #include <urkbias.h>
61
GatherCDSNew(void)62 extern Gather_CDSPtr GatherCDSNew (void)
63 {
64 Gather_CDSPtr gcdsp;
65
66 if ((gcdsp = (Gather_CDSPtr) MemNew (sizeof (Gather_CDS))) == NULL)
67 return NULL;
68 gcdsp->LOscore = -1.0;
69 gcdsp->scorecut = 0.5;
70
71 return gcdsp;
72 }
73
GatherCDSFree(Gather_CDSPtr gcdsp)74 extern Gather_CDSPtr GatherCDSFree (Gather_CDSPtr gcdsp)
75 {
76 SeqLocPtr slp, slpn;
77 SeqIdPtr id;
78
79 if (gcdsp == NULL)
80 return NULL;
81
82 MemFree (gcdsp->tableGlobal);
83 MemFree (gcdsp->tableRefine);
84
85 slp = gcdsp->slpGlobal;
86 while (slp != NULL)
87 {
88 slpn = slp->next;
89 id = SeqLocId (slp);
90 if (id != NULL)
91 id->next = SeqIdSetFree (id->next);
92 SeqLocFree (slp);
93 slp = slpn;
94 }
95 slp = gcdsp->slpRefine;
96 while (slp != NULL)
97 {
98 slpn = slp->next;
99 id = SeqLocId (slp);
100 if (id != NULL)
101 id->next = SeqIdSetFree (id->next);
102 SeqLocFree (slp);
103 slp = slpn;
104 }
105 slp = gcdsp->slpAll;
106 while (slp != NULL)
107 {
108 slpn = slp->next;
109 id = SeqLocId (slp);
110 if (id != NULL)
111 id->next = SeqIdSetFree (id->next);
112 SeqLocFree (slp);
113 slp = slpn;
114 }
115 slp = gcdsp->slpHit;
116 while (slp != NULL)
117 {
118 slpn = slp->next;
119 id = SeqLocId (slp);
120 if (id != NULL)
121 id->next = SeqIdSetFree (id->next);
122 SeqLocFree (slp);
123 slp = slpn;
124 }
125 slp = gcdsp->slpFound;
126 while (slp != NULL)
127 {
128 slpn = slp->next;
129 id = SeqLocId (slp);
130 if (id != NULL)
131 id->next = SeqIdSetFree (id->next);
132 SeqLocFree (slp);
133 slp = slpn;
134 }
135
136 return (Gather_CDSPtr) MemFree (gcdsp);
137 }
138
SeqLocDup(SeqLocPtr slpold)139 extern SeqLocPtr SeqLocDup (SeqLocPtr slpold)
140 {
141 SeqLocPtr slpnew, slpn, slp;
142 SeqIdPtr id;
143
144 if (slpold == NULL)
145 return NULL;
146
147 slpnew = (SeqLocPtr) AsnIoMemCopy ((Pointer) slpold,
148 (AsnReadFunc) SeqLocAsnRead,
149 (AsnWriteFunc) SeqLocAsnWrite);
150 slp = slpnew->next;
151 while (slp != NULL)
152 {
153 slpn = slp->next;
154 id = SeqLocId (slp);
155 if (id != NULL)
156 id->next = SeqIdSetFree (id->next);
157 SeqLocFree (slp);
158 slp = slpn;
159 }
160 slpnew->next = NULL;
161
162 return slpnew;
163 }
164
SeqLocDupAll(SeqLocPtr slpold)165 extern SeqLocPtr SeqLocDupAll (SeqLocPtr slpold)
166 {
167 SeqLocPtr slpnew;
168
169 if (slpold == NULL)
170 return NULL;
171
172 slpnew = (SeqLocPtr) AsnIoMemCopy ((Pointer) slpold,
173 (AsnReadFunc) SeqLocAsnRead,
174 (AsnWriteFunc) SeqLocAsnWrite);
175 return slpnew;
176 }
177
SeqLocLink(SeqLocPtr PNTR slph,SeqLocPtr slp)178 extern SeqLocPtr SeqLocLink (SeqLocPtr PNTR slph, SeqLocPtr slp)
179 {
180 SeqLocPtr slpnext;
181
182 if (slph == NULL || *slph == NULL)
183 {
184 *slph = slp;
185 return *slph;
186 }
187 slpnext = *slph;
188 while (slpnext->next != NULL)
189 slpnext = slpnext->next;
190 slpnext->next = slp;
191 return *slph;
192 }
193
SeqLocMatch(SeqLocPtr slplist,SeqLocPtr slpnew)194 extern Boolean SeqLocMatch (SeqLocPtr slplist, SeqLocPtr slpnew)
195 {
196 Int4 start1, stop1, start2, stop2;
197
198 if (slpnew == NULL || slplist == NULL)
199 {
200 return FALSE;
201 }
202
203 while (slplist != NULL)
204 {
205 start1 = SeqLocStart (slpnew);
206 stop1 = SeqLocStop (slpnew);
207 start2 = SeqLocStart (slplist);
208 stop2 = SeqLocStop (slplist);
209
210 if (start1 == start2 && stop1 == stop2)
211 return TRUE;
212 slplist = slplist->next;
213 }
214 return FALSE;
215 }
216
UniqueOrfs(SeqLocPtr PNTR pslpFound)217 extern void UniqueOrfs (SeqLocPtr PNTR pslpFound)
218 {
219 SeqLocPtr slpFound, slpNew, slp;
220 SeqIdPtr id;
221
222 slpFound = *pslpFound;
223 slpNew = NULL;
224 while (slpFound != NULL)
225 {
226 if (!SeqLocMatch (slpNew, slpFound))
227 SeqLocLink (&slpNew, SeqLocDup (slpFound));
228 slpFound = slpFound->next;
229 }
230 slpFound = *pslpFound;
231 while (slpFound != NULL)
232 {
233 slp = slpFound->next;
234 id = SeqLocId (slpFound);
235 if (id != NULL)
236 id->next = SeqIdSetFree (id->next);
237 SeqLocFree (slpFound);
238 slpFound = slp;
239 }
240 *pslpFound = slpNew;
241 return;
242 }
243
SeqLocCompProc(VoidPtr ptr1,VoidPtr ptr2)244 static int LIBCALLBACK SeqLocCompProc (VoidPtr ptr1, VoidPtr ptr2)
245 {
246 SeqLocPtr slp1, slp2;
247
248 if (ptr1 != NULL && ptr2 != NULL)
249 {
250 slp1 = *((SeqLocPtr PNTR) ptr1);
251 slp2 = *((SeqLocPtr PNTR) ptr2);
252 if (slp1 != NULL && slp2 != NULL)
253 {
254 if (SeqLocStart (slp1) > SeqLocStart (slp2))
255 {
256 return 1;
257 }
258 else if (SeqLocStart (slp1) < SeqLocStart (slp2))
259 {
260 return -1;
261 }
262 }
263 }
264 return 0;
265 }
266
SortByLocOrfs(SeqLocPtr PNTR pslpFound)267 extern void SortByLocOrfs (SeqLocPtr PNTR pslpFound)
268 {
269 SeqLocPtr slpFound, PNTR slpt;
270 Int4 num, i;
271
272 num = 0;
273 slpFound = *pslpFound;
274 while (slpFound != NULL)
275 {
276 num++;
277 slpFound = slpFound->next;
278 }
279
280 slpt = MemNew ((size_t)(num+1) * sizeof (SeqLocPtr));
281 slpFound = *pslpFound;
282 i = 0;
283 while (slpFound != NULL)
284 {
285 slpt[i] = slpFound;
286 slpFound = slpFound->next;
287 i++;
288 }
289
290 HeapSort (slpt, num, sizeof (SeqLocPtr), SeqLocCompProc);
291
292 for (i = 0; i < num; i++)
293 {
294 slpFound = slpt[i];
295 slpFound->next = slpt[i+1];
296 }
297
298 *pslpFound = slpt[0];
299 MemFree (slpt);
300
301 return;
302 }
303
RemoveInternalOrfs(SeqLocPtr PNTR slpFound)304 extern void RemoveInternalOrfs (SeqLocPtr PNTR slpFound)
305 {
306 SeqLocPtr slp1, slp2, slp = NULL;
307 SeqIdPtr id;
308 Int4 start1, stop1, start2, stop2;
309 Boolean flagInternal;
310
311 slp1 = *slpFound;
312 while (slp1 != NULL)
313 {
314 start1 = SeqLocStart (slp1);
315 stop1 = SeqLocStop (slp1);
316 flagInternal = FALSE;
317 slp2 = *slpFound;
318 while (slp2 != NULL)
319 {
320 start2 = SeqLocStart (slp2);
321 stop2 = SeqLocStop (slp2);
322 if ((start1 > start2 && start1 < stop2) &&
323 (stop1 > start2 && stop1 < stop2))
324 {
325 flagInternal = TRUE;
326 break;
327 }
328 slp2 = slp2->next;
329 }
330 if (!flagInternal)
331 SeqLocLink (&slp, SeqLocDup (slp1));
332 slp1 = slp1->next;
333 }
334 slp1 = *slpFound;
335 while (slp1 != NULL)
336 {
337 slp2 = slp1->next;
338 id = SeqLocId (slp1);
339 if (id != NULL)
340 id->next = SeqIdSetFree (id->next);
341 SeqLocFree (slp1);
342 slp1 = slp2;
343 }
344 *slpFound = slp;
345 return;
346 }
347
ScanCodonUsage(GatherContextPtr gcp)348 static Boolean ScanCodonUsage (GatherContextPtr gcp)
349 {
350 Gather_CDSPtr gcdsp;
351 BioseqPtr bsp;
352 SeqLocPtr slp, slpt;
353 Int4Ptr cutp;
354 FloatHi score;
355
356 if (gcp == NULL)
357 return FALSE;
358 if ((gcdsp = (Gather_CDSPtr) gcp->userdata) == NULL)
359 return FALSE;
360 if (gcp->thistype != OBJ_BIOSEQ)
361 return TRUE;
362 if ((bsp = (BioseqPtr) (gcp->thisitem)) == NULL)
363 return TRUE;
364
365 if (gcdsp->findcount != 0)
366 return TRUE;
367
368 slp = gcdsp->slpAll;
369 while (slp != NULL)
370 {
371 cutp = NewCodonTable ();
372 slpt = SeqLocDup (slp);
373 AddSeqLocToCodonTable (cutp, bsp, slpt, TRUE);
374 score = Confide (cutp, gcdsp->tableRefine);
375 FreeCodonTable (cutp);
376 if (score < gcdsp->scorecut)
377 {
378 SeqLocLink (&(gcdsp->slpHit), SeqLocDup (slp));
379 gcdsp->findcount++;
380 }
381 SeqLocFree (slpt);
382 slp = slp->next;
383 }
384 return TRUE;
385 }
386
RefineCodonUsage(GatherContextPtr gcp)387 static Boolean RefineCodonUsage (GatherContextPtr gcp)
388 {
389 Gather_CDSPtr gcdsp;
390 BioseqPtr bsp;
391 SeqLocPtr slp, slpt;
392 Int4Ptr cutp;
393 FloatHi score;
394
395 if (gcp == NULL)
396 return FALSE;
397 if ((gcdsp = (Gather_CDSPtr) gcp->userdata) == NULL)
398 return FALSE;
399 if (gcp->thistype != OBJ_BIOSEQ)
400 return TRUE;
401 if ((bsp = (BioseqPtr) (gcp->thisitem)) == NULL)
402 return TRUE;
403
404 if (gcdsp->tableRefine != NULL)
405 return TRUE;
406
407 slp = gcdsp->slpAll;
408 while (slp != NULL)
409 {
410 cutp = NewCodonTable ();
411 slpt = SeqLocDup (slp);
412 AddSeqLocToCodonTable (cutp, bsp, slpt, TRUE);
413 score = Confide (cutp, gcdsp->tableGlobal);
414 FreeCodonTable (cutp);
415 if (score < gcdsp->scorecut)
416 {
417 gcdsp->refinecount++;
418 if (gcdsp->tableRefine == NULL)
419 gcdsp->tableRefine = NewCodonTable ();
420 AddSeqLocToCodonTable (gcdsp->tableRefine, bsp, slpt, TRUE);
421 SeqLocLink (&(gcdsp->slpRefine), SeqLocDup (slp));
422 }
423 SeqLocFree (slpt);
424 slp = slp->next;
425 }
426 return TRUE;
427 }
428
StandardMean(GatherContextPtr gcp)429 static Boolean StandardMean (GatherContextPtr gcp)
430 {
431 Gather_CDSPtr gcdsp;
432 BioseqPtr bsp;
433 SeqLocPtr slp, slpt;
434 Int4Ptr cutp;
435 FloatHi score;
436
437 if (gcp == NULL)
438 return FALSE;
439 if ((gcdsp = (Gather_CDSPtr) gcp->userdata) == NULL)
440 return FALSE;
441 if (gcp->thistype != OBJ_BIOSEQ)
442 return TRUE;
443 if ((bsp = (BioseqPtr) (gcp->thisitem)) == NULL)
444 return TRUE;
445
446 if (gcdsp->stdcount != 0)
447 return TRUE;
448
449 slp = gcdsp->slpAll;
450 while (slp != NULL)
451 {
452 cutp = NewCodonTable ();
453 slpt = SeqLocDup (slp);
454 AddSeqLocToCodonTable (cutp, bsp, slpt, TRUE);
455 score = Confide (cutp, gcdsp->tableGlobal);
456 FreeCodonTable (cutp);
457 gcdsp->mean += score;
458 gcdsp->stdcount++;
459 if (gcdsp->LOscore == -1)
460 gcdsp->LOscore = score;
461 if (score < gcdsp->LOscore)
462 gcdsp->LOscore = score;
463 if (score > gcdsp->HIscore)
464 gcdsp->HIscore = score;
465 SeqLocFree (slpt);
466 slp = slp->next;
467 }
468 return TRUE;
469 }
470
StandardDeviation(GatherContextPtr gcp)471 static Boolean StandardDeviation (GatherContextPtr gcp)
472 {
473 Gather_CDSPtr gcdsp;
474 BioseqPtr bsp;
475 SeqLocPtr slp, slpt;
476 Int4Ptr cutp;
477 FloatHi score;
478
479 if (gcp == NULL)
480 return FALSE;
481 if ((gcdsp = (Gather_CDSPtr) gcp->userdata) == NULL)
482 return FALSE;
483 if (gcp->thistype != OBJ_BIOSEQ)
484 return TRUE;
485 if ((bsp = (BioseqPtr) (gcp->thisitem)) == NULL)
486 return TRUE;
487
488 if (gcdsp->stdev != 0.0)
489 return TRUE;
490
491 slp = gcdsp->slpAll;
492 while (slp != NULL)
493 {
494 cutp = NewCodonTable ();
495 slpt = SeqLocDup (slp);
496 AddSeqLocToCodonTable (cutp, bsp, slpt, TRUE);
497 score = Confide (cutp, gcdsp->tableGlobal);
498 FreeCodonTable (cutp);
499 score -= gcdsp->mean;
500 score *= score;
501 gcdsp->stdev += score;
502 SeqLocFree (slpt);
503 slp = slp->next;
504 }
505 return TRUE;
506 }
507
MakeGlobalTable(GatherContextPtr gcp)508 static Boolean MakeGlobalTable (GatherContextPtr gcp)
509 {
510 Gather_CDSPtr gcdsp;
511 BioseqPtr bsp;
512 SeqLocPtr slp, slpt;
513
514 if (gcp == NULL)
515 return FALSE;
516 if ((gcdsp = (Gather_CDSPtr) gcp->userdata) == NULL)
517 return FALSE;
518 if (gcp->thistype != OBJ_BIOSEQ)
519 return TRUE;
520 if ((bsp = (BioseqPtr) (gcp->thisitem)) == NULL)
521 return TRUE;
522
523 if (gcdsp->tableGlobal != NULL)
524 return TRUE;
525
526 slp = gcdsp->slpGlobal;
527 while (slp != NULL)
528 {
529 gcdsp->globalcount++;
530 if (gcdsp->tableGlobal == NULL)
531 gcdsp->tableGlobal = NewCodonTable ();
532 slpt = SeqLocDup (slp);
533 AddSeqLocToCodonTable (gcdsp->tableGlobal, bsp, slpt, TRUE);
534 SeqLocFree (slpt);
535 slp = slp->next;
536 }
537 return TRUE;
538 }
539
CullGlobalOrfs(SeqLocPtr PNTR pslpGlobal,SeqLocPtr PNTR pslpRefine)540 static Int4 CullGlobalOrfs (SeqLocPtr PNTR pslpGlobal,
541 SeqLocPtr PNTR pslpRefine)
542 {
543 Int4 globalcount;
544 SeqLocPtr slpGlobal, slpRefine, slpNew, slp;
545
546 slpGlobal = *pslpGlobal;
547 slpRefine = *pslpRefine;
548 slpNew = NULL;
549 globalcount = 0;
550 while (slpGlobal != NULL)
551 {
552 if (!SeqLocMatch (slpRefine, slpGlobal))
553 {
554 globalcount++;
555 SeqLocLink (&slpNew, SeqLocDup (slpGlobal));
556 }
557 slpGlobal = slpGlobal->next;
558 }
559 slpGlobal = *pslpGlobal;
560 while (slpGlobal != NULL)
561 {
562 slp = slpGlobal->next;
563 SeqLocFree (slpGlobal);
564 slpGlobal = slp;
565 }
566 *pslpGlobal = slpNew;
567 return globalcount;
568 }
569
FindSimilarBiasOrfs(SeqEntryPtr sep,FloatHi probcut,Int4 clustmin,Int4 findmin,SeqLocPtr slpKnown,SeqLocPtr slpPotential)570 extern SeqLocPtr FindSimilarBiasOrfs (SeqEntryPtr sep, FloatHi probcut,
571 Int4 clustmin, Int4 findmin,
572 SeqLocPtr slpKnown,
573 SeqLocPtr slpPotential)
574 {
575 static GatherScope gs;
576 GatherScopePtr gsp;
577 Gather_CDSPtr gcdsp;
578
579 Int4 gcount;
580 SeqLocPtr slp, slpn;
581
582 if (probcut == 0.0)
583 probcut = 0.5;
584 if (clustmin == 0)
585 clustmin = 2;
586 if (findmin == 0)
587 findmin = 4;
588
589 gsp = &gs;
590 gcdsp = GatherCDSNew ();
591
592 MemSet ((Pointer) gsp, 0, sizeof (GatherScope));
593 MemSet ((Pointer) gsp->ignore, (int) (TRUE),
594 (size_t) (OBJ_MAX * sizeof (Boolean)));
595 gsp->ignore[OBJ_BIOSEQ] = FALSE;
596
597 slp = slpKnown;
598 while (slp != NULL)
599 {
600 SeqLocLink (&gcdsp->slpGlobal, SeqLocDup (slp));
601 slp = slp->next;
602 }
603 slp = slpPotential;
604 while (slp != NULL)
605 {
606 SeqLocLink (&gcdsp->slpAll, SeqLocDup (slp));
607 slp = slp->next;
608 }
609
610 GatherSeqEntry (sep, (Pointer) gcdsp, MakeGlobalTable, (Pointer) gsp);
611
612 while (gcdsp->tableGlobal != NULL)
613 {
614 gcdsp->stdcount = 0;
615 gcdsp->refinecount = 0;
616 gcdsp->findcount = 0;
617 gcdsp->stdev = 0.0;
618 gcdsp->HIscore = 0.0;
619 gcdsp->LOscore = -1.0;
620 gcdsp->mean = 0.0;
621 gcdsp->scorecut = 0.5;
622
623 GatherSeqEntry (sep, (Pointer) gcdsp, StandardMean, (Pointer) gsp);
624 if (gcdsp->stdcount > 0)
625 gcdsp->mean /= gcdsp->stdcount;
626 GatherSeqEntry (sep, (Pointer) gcdsp, StandardDeviation, (Pointer) gsp);
627
628 if (gcdsp->stdcount > 1)
629 gcdsp->stdev /= (gcdsp->stdcount - 1);
630 else
631 gcdsp->stdev = 0.0;
632 gcdsp->stdev = (FloatHi) sqrt (gcdsp->stdev);
633 gcdsp->scorecut = gcdsp->LOscore + (gcdsp->stdev * probcut);
634
635 slp = gcdsp->slpRefine;
636 while (slp != NULL)
637 {
638 slpn = slp->next;
639 SeqLocFree (slp);
640 slp = slpn;
641 }
642 gcdsp->slpRefine = slp;
643
644 gcdsp->tableRefine = FreeCodonTable (gcdsp->tableRefine);
645 GatherSeqEntry (sep, (Pointer) gcdsp, RefineCodonUsage, (Pointer) gsp);
646 if (gcdsp->tableRefine != NULL)
647 {
648 if (gcdsp->refinecount >= clustmin)
649 {
650 gcdsp->scorecut *= 1.5; /* increase a bit to see any branch jumps */
651 GatherSeqEntry (sep, (Pointer) gcdsp, ScanCodonUsage, (Pointer) gsp);
652 }
653 if (gcdsp->findcount < findmin)
654 {
655 slp = gcdsp->slpHit;
656 while (slp != NULL)
657 {
658 slpn = slp->next;
659 SeqLocFree (slp);
660 slp = slpn;
661 }
662 gcdsp->slpHit = slp;
663 }
664 else
665 {
666 SeqLocLink (&gcdsp->slpFound, gcdsp->slpHit);
667 gcdsp->slpHit = NULL;
668 }
669 gcount = CullGlobalOrfs (&gcdsp->slpGlobal, &gcdsp->slpRefine);
670 gcdsp->tableGlobal = FreeCodonTable (gcdsp->tableGlobal);
671 if (gcount != gcdsp->globalcount)
672 {
673 gcdsp->globalcount = 0;
674 GatherSeqEntry (sep, (Pointer) gcdsp, MakeGlobalTable, (Pointer) gsp);
675 }
676 }
677 else
678 {
679 slp = gcdsp->slpGlobal;
680 while (slp != NULL)
681 {
682 slpn = slp->next;
683 SeqLocFree (slp);
684 slp = slpn;
685 }
686 gcdsp->slpGlobal = slp;
687 gcdsp->tableGlobal = FreeCodonTable (gcdsp->tableGlobal);
688 }
689 }
690 UniqueOrfs (&gcdsp->slpFound);
691 slp = gcdsp->slpFound;
692 gcdsp->slpFound = NULL;
693 GatherCDSFree (gcdsp);
694 return slp;
695 }
696
BiasScoreBioseq(BioseqPtr bsp,Int4Ptr tableGlobal,Int4 tripletwindow,Int4 xframe,Uint1 xstrand)697 extern FloatHiPtr BiasScoreBioseq (BioseqPtr bsp, Int4Ptr tableGlobal,
698 Int4 tripletwindow, Int4 xframe,
699 Uint1 xstrand)
700 {
701 FloatHiPtr score;
702 Int4 iscore;
703
704 SeqLocPtr slp;
705 SeqIntPtr sint;
706 Int4 start, stop, xstop, xwindow;
707 Int4Ptr cutp;
708
709 if (bsp == NULL)
710 return NULL;
711 if (!ISA_na (bsp->mol))
712 return NULL;
713 if (bsp->length < tripletwindow)
714 return NULL;
715
716 slp = ValNodeNew (NULL);
717 sint = SeqIntNew ();
718 slp->choice = SEQLOC_INT;
719 slp->data.ptrvalue = sint;
720
721 xwindow = tripletwindow;
722 xstop = (bsp->length + 3 - xframe - xwindow) / 3;
723 score = (FloatHiPtr) MemNew ((size_t) (sizeof (FloatHi) * xstop));
724 xwindow -= 3;
725 xstop--;
726
727 start = xframe;
728 stop = start + xwindow - 1;
729 sint->from = start;
730 sint->to = stop;
731 sint->strand = xstrand;
732 cutp = CodonTableFromSeqLoc (bsp, slp);
733
734 iscore = 0;
735
736 xstop = bsp->length - 3;
737 for (start = stop + 1; start <= xstop; start += 3)
738 {
739 sint->from = start;
740 sint->to = start + 2;
741 AddSeqLocToCodonTable (cutp, bsp, slp, TRUE);
742 score[iscore++] = Confide (cutp, tableGlobal);
743 sint->from -= xwindow;
744 sint->to -= xwindow;
745 AddSeqLocToCodonTable (cutp, bsp, slp, FALSE);
746 }
747
748 FreeCodonTable (cutp);
749 slp->data.ptrvalue = (Pointer) SeqIntFree (sint);
750 SeqLocFree (slp);
751 return score;
752 }
753