← Back to branch summary

~ubuntu-branches/ubuntu/breezy/ncbi-tools6/breezy

« back to all changes in this revision

Viewing changes to tools/urkbias.c

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

Show diffs side-by-side

added added

removed removed

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