1 #ifdef _cplusplus
2 extern "C" {
3 #endif
4 #include "aligngenemodel.h"
5
6
7
8
9 # line 77 "aligngenemodel.dy"
weighted_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)10 Probability weighted_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)
11 {
12 int aa;
13 int seq_no;
14 Probability out;
15 Probability aa_prob;
16
17 int codon;
18
19
20 out = 0.0;
21
22 for(aa=0;aa<26;aa++) {
23 if( agmp->rm->aminoacid[aa] < 0.0000000001 ) {
24 continue;
25 }
26
27 aa_prob = agmp->rm->aminoacid[aa];
28
29 for(seq_no=0;seq_no<store->len;seq_no++) {
30
31 if( !isalpha(store->codon[seq_no].seq[0]) ||
32 !isalpha(store->codon[seq_no].seq[1]) ||
33 !isalpha(store->codon[seq_no].seq[2]) ) {
34 return 0.5;
35 }
36
37 codon = codon_from_seq(store->codon[seq_no].seq);
38 if( is_stop_codon(codon,agmp->ct) ) {
39 if( seq_no > 0 && agmp->tolerate_nonanchor_stops ) {
40 return agmp->nonanchor_stop;
41 } else {
42 return 0.0;
43 }
44 }
45
46 aa_prob *= agmp->protein->comp[aa][aminoacid_from_codon(agmp->ct,codon)-'A'] * store->codon[seq_no].weight;
47 }
48
49 out += aa_prob;
50 }
51
52 return out;
53 }
54
55 # line 122 "aligngenemodel.dy"
window_AlignGeneModel(SeqAlign * sal,AlignGeneModel * agm,AlignGeneModelParam * agmp)56 void window_AlignGeneModel(SeqAlign * sal,AlignGeneModel * agm,AlignGeneModelParam * agmp)
57 {
58 int i;
59 int j;
60 int k;
61 int count;
62 int win_count;
63 double average_change;
64 int change;
65
66 double * change_array;
67
68 change_array = calloc(agm->len,sizeof(double));
69
70
71 fprintf(stderr,"Got Window at %f and non window at %f\n",agmp->coding_window_bonus,agmp->noncoding_window_pen);
72
73 for(i=1;i<sal->len;i++) {
74
75 average_change = 0.0;
76 count = 0;
77
78 for(j=0;j<sal->seq[0]->len-32;j++) {
79 change = 0;
80 win_count = 0;
81 for(k = 0;k<31;k++) {
82 if( sal->seq[0]->seq[j+k] != sal->seq[i]->seq[j+k] && sal->seq[i]->seq[j+k] != '~') {
83 change++;
84 }
85 if( sal->seq[i]->seq[j+k] != '~' ) {
86 win_count++;
87 }
88 }
89 if( win_count == 0 ) {
90 change_array[j+15] = 0.9;
91 continue;
92 }
93
94 change_array[j+15] = ((double) change) / win_count;
95 average_change += change_array[j+15];
96 count++;
97 /* fprintf(stdout,"Got %f, running average %f\n",change_array[j+15],average_change/(double)count);*/
98 }
99
100 average_change = average_change / (double) count;
101
102 fprintf(stdout,"For %s, Average change is %f\n",sal->seq[i]->name,average_change);
103
104 for(j=15;j<sal->seq[0]->len-16;j++) {
105
106 if( change_array[j] < (average_change/agmp->coding_window_thres) ) {
107 /* fprintf(stdout,"At position %d in %s giving coding bonus\n",j,sal->seq[i]->name);*/
108
109 agm->forward_coding[j] *= agmp->coding_window_bonus;
110 agm->reverse_coding[j] *= agmp->coding_window_bonus;
111 }
112 if( change_array[j] > (average_change * agmp->noncoding_window_thres) ) {
113 /* fprintf(stdout,"At position %d in %s giving penalty\n",j,sal->seq[i]->name);*/
114 agm->forward_coding[j] *= agmp->noncoding_window_pen;
115 agm->reverse_coding[j] *= agmp->noncoding_window_pen;
116 }
117 }
118 }
119
120 free(change_array);
121
122 return;
123 }
124
125
126 # line 192 "aligngenemodel.dy"
create_NonCodingSimpleModel(Probability change)127 NonCodingSimpleModel * create_NonCodingSimpleModel(Probability change)
128 {
129 NonCodingSimpleModel * out;
130
131 out = NonCodingSimpleModel_alloc();
132
133 out->one_off = change * (1.0-change) * (1.0-change) * 3;
134 out->two_off = change * change * (1.0-change) * 3;
135 out->identical = 1.0 - (out->one_off + out->two_off);
136
137 return out;
138 }
139
140 # line 205 "aligngenemodel.dy"
simple_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)141 Probability simple_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)
142 {
143 int i;
144 int j;
145 int k;
146 Probability out = 1.0;
147 Probability col = 0.0;
148 Probability c = 0.0;
149
150
151 for(k=0;k<3;k++) {
152 col = 0.0;
153 for(j=0;j<4;j++) {
154 c= 0.25;
155 for(i=0;i<store->len;i++) {
156 if( store->codon[i].seq[2-k] == '-' || store->codon[i].seq[2-k] == '~') {
157 return 0.5;
158 }
159 c *= agmp->dm->prob[j][base_from_char(store->codon[i].seq[2-k])];
160 }
161 col += c;
162 }
163
164 out *= col;
165 }
166
167 return out;
168 }
169
170 # line 234 "aligngenemodel.dy"
weighted_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)171 Probability weighted_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)
172 {
173 int i;
174 int j;
175 int k;
176
177
178 Probability out = 1.0;
179 Probability col = 0.0;
180 Probability c = 0.0;
181
182
183 for(k=0;k<3;k++) {
184 col = 0.0;
185 for(j=0;j<4;j++) {
186 c= 0.25;
187 for(i=0;i<store->len;i++) {
188 if( store->codon[i].seq[2-k] == '-' || store->codon[i].seq[2-k] == '~') {
189 return 0.5;
190 }
191 c *= agmp->dm->prob[j][base_from_char(store->codon[i].seq[2-k])];
192 }
193 col += c;
194 }
195
196 out *= col;
197 }
198
199 return out;
200
201
202 }
203
204
205 # line 268 "aligngenemodel.dy"
weighted_simple_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)206 Probability weighted_simple_non_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)
207 {
208 int i;
209 int j;
210 int k;
211
212
213 Probability out = 1.0;
214 Probability col = 0.0;
215 Probability c = 0.0;
216
217
218 for(k=0;k<3;k++) {
219 col = 0.0;
220 for(j=0;j<4;j++) {
221 c= 0.25;
222 for(i=0;i<store->len;i++) {
223 if( store->codon[i].seq[2-k] == '-' ) {
224 return 0.5;
225 }
226 c *= agmp->dm->prob[j][base_from_char(store->codon[i].seq[2-k])];
227 }
228 col += c;
229 }
230
231 out *= col;
232 }
233
234 return out;
235
236
237 }
238
239 # line 301 "aligngenemodel.dy"
simple_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)240 Probability simple_coding_AlignGeneColumn(AlignGeneColumnStore * store,AlignGeneModelParam * agmp)
241 {
242 int i;
243 int j;
244 int codon;
245 Probability out = 0.0;
246 Probability c = 0.0;
247 Probability f;
248
249 assert(store);
250 assert(agmp);
251
252 for(j=0;j<26;j++) {
253 if( agmp->rm->aminoacid[j] < 0.0000000001 ) {
254 continue;
255 }
256 c = agmp->rm->aminoacid[j];
257 for(i=0;i<store->len;i++) {
258 if( !isalpha(store->codon[i].seq[0]) ||
259 !isalpha(store->codon[i].seq[1]) ||
260 !isalpha(store->codon[i].seq[2]) ) {
261 return 0.0;
262 }
263
264 codon = codon_from_seq(store->codon[i].seq);
265 if( is_stop_codon(codon,agmp->ct) ) {
266 return 0.0;
267 }
268 f = agmp->protein->comp[j][aminoacid_from_codon(agmp->ct,codon)-'A'];
269
270 /*
271 fprintf(stdout,"For amino acid %c, match %.4f cum %.4f\n",j+'A',
272 agmp->protein->comp[j][aminoacid_from_codon(agmp->ct,codon)-'A'],
273 c);
274 */
275
276 /* c = c*f;*/
277
278 c = c*f;
279 /* fprintf(stdout,"Refactoring by %f\n",store->codon[i].weight); */
280 /* c = c*100; */
281 }
282 out += c;
283 }
284
285 /* fprintf(stdout,"final %.4f\n",out);*/
286 return out;
287 }
288
289 # line 350 "aligngenemodel.dy"
fill_reverse_AlignGeneColumn(AlignGeneColumnStore * store,SeqAlign * sa,int forward_codon_end)290 boolean fill_reverse_AlignGeneColumn(AlignGeneColumnStore * store,SeqAlign * sa,int forward_codon_end)
291 {
292 int i;
293
294 assert(store);
295 assert(sa);
296 assert(store->len == sa->len);
297 assert(forward_codon_end >= 0 && forward_codon_end < sa->seq[0]->len);
298
299 for(i=0;i<sa->len;i++) {
300 store->codon[i].seq[0] = char_complement_base(sa->seq[i]->seq[forward_codon_end]);
301 store->codon[i].seq[1] = char_complement_base(sa->seq[i]->seq[forward_codon_end-1]);
302 store->codon[i].seq[2] = char_complement_base(sa->seq[i]->seq[forward_codon_end-2]);
303 store->codon[i].weight = sa->seq[i]->weight;
304 }
305
306 return TRUE;
307
308 }
309
310 # line 370 "aligngenemodel.dy"
fill_forward_AlignGeneColumn(AlignGeneColumnStore * store,SeqAlign * sa,int codon_end)311 boolean fill_forward_AlignGeneColumn(AlignGeneColumnStore * store,SeqAlign * sa,int codon_end)
312 {
313 int i;
314
315 assert(store);
316 assert(sa);
317 assert(store->len == sa->len);
318 assert(codon_end >= 0 && codon_end < sa->seq[0]->len);
319
320 for(i=0;i<sa->len;i++) {
321 store->codon[i].seq[0] = sa->seq[i]->seq[codon_end-2];
322 store->codon[i].seq[1] = sa->seq[i]->seq[codon_end-1];
323 store->codon[i].seq[2] = sa->seq[i]->seq[codon_end];
324 store->codon[i].weight = sa->seq[i]->weight;
325 }
326
327 return TRUE;
328 }
329
330 # line 389 "aligngenemodel.dy"
new_empty_AlignGeneColumnStore(int len)331 AlignGeneColumnStore * new_empty_AlignGeneColumnStore(int len)
332 {
333 AlignGeneColumnStore * out;
334
335 out = AlignGeneColumnStore_alloc();
336
337 out->codon = calloc(len,sizeof(AlignGeneCodon));
338 out->len = len;
339
340 return out;
341 }
342
343 # line 401 "aligngenemodel.dy"
free_AlignGeneCodon(AlignGeneCodon * c)344 AlignGeneCodon * free_AlignGeneCodon(AlignGeneCodon * c)
345 {
346 ckfree(c);
347 return NULL;
348 }
349
350
351 # line 408 "aligngenemodel.dy"
new_AlignGeneModel(int len)352 AlignGeneModel * new_AlignGeneModel(int len)
353 {
354 AlignGeneModel * out;
355
356 out = AlignGeneModel_alloc();
357
358 out->len = len;
359
360 out->forward_coding = calloc(len,sizeof(Probability));
361 out->reverse_coding = calloc(len,sizeof(Probability));
362 out->splice5_forward = calloc(len,sizeof(Probability));
363 out->splice3_forward = calloc(len,sizeof(Probability));
364 out->splice5_reverse = calloc(len,sizeof(Probability));
365 out->splice3_reverse = calloc(len,sizeof(Probability));
366 out->change_rate = calloc(len,sizeof(double));
367
368 out->align = NULL;
369 out->anchor = NULL;
370
371 return out;
372 }
373
374 # line 430 "aligngenemodel.dy"
show_AlignGeneModel(AlignGeneModel * agm,SeqAlign * sal,CodonTable * ct,GenomicRegion * gr,FILE * ofp,AlignGeneModelParam * agmp)375 void show_AlignGeneModel(AlignGeneModel * agm,SeqAlign * sal,CodonTable * ct,GenomicRegion * gr,FILE * ofp,AlignGeneModelParam * agmp)
376 {
377 int i;
378 int j;
379 int * should_show;
380 int ex_len;
381 int k;
382 AlignGeneColumnStore * store;
383 int tr_len;
384 int phase_score[3];
385 int total_phase_0 = 0;
386 int phase;
387 int * phase_array;
388 int pos;
389 int count = 0;
390
391 fprintf(stderr,"Gene model at %d with %d\n",gr,gr->len);
392
393 should_show = calloc(agm->len,sizeof(int));
394 phase_array = calloc(agm->len,sizeof(int));
395
396 fprintf(stderr,"Gene model at %d with %d\n",gr,gr->len);
397
398 if( gr != NULL ) {
399 if( gr->len == 0 ) {
400 warn("genomic region given, but no genes on region!");
401 }
402
403 for(i=0;i<agm->len;i++)
404 should_show[i] = 0;
405 } else {
406 for(i=0;i<agm->len;i++)
407 should_show[i] = 1;
408 }
409
410
411
412 if( gr != NULL ) {
413 for( i=0;i<gr->len;i++) {
414 tr_len = 3;
415 for( j=0;j<gr->gene[i]->transcript[0]->ex_len; j++) {
416 ex_len = gr->gene[i]->transcript[0]->exon[j]->end -
417 gr->gene[i]->transcript[0]->exon[j]->start;
418
419 phase_score[0] = 0;
420 phase_score[1] = 0;
421 phase_score[2] = 0;
422
423
424 /* deal with splice sites etc */
425 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start-4;
426 should_show[pos] = 1;
427 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start-3;
428 should_show[pos] = 1;
429 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start-2;
430 should_show[pos] = 1;
431 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start-1;
432 should_show[pos] = 1;
433
434 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start+ex_len+1;
435 should_show[pos] = 1;
436 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start+ex_len;
437 should_show[pos] = 1;
438 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start+ex_len+2;
439 should_show[pos] = 1;
440 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start+ex_len+3;
441 should_show[pos] = 1;
442
443 for(k=0;k<ex_len;k++) {
444 phase = ((tr_len -2) %3);
445 pos = gr->gene[i]->start+gr->gene[i]->transcript[0]->exon[j]->start+k;
446
447 phase_score[phase] += Probability2Score(agm->forward_coding[pos]);
448 phase_array[pos] = phase_score[phase];
449 /*
450 printf("%d Adding %d in to %d with %f\n",phase,Probability2Score(agm->forward_coding[pos]),phase_score[phase],Score2Bits(phase_score[phase]));
451 */
452
453 if( phase == 0 ) {
454 total_phase_0 += Probability2Score(agm->forward_coding[i]);
455 should_show[pos] = 2;
456 } else {
457 should_show[pos] = 1;
458 }
459 tr_len++;
460 }
461 }
462 }
463 }
464
465 store = new_empty_AlignGeneColumnStore(sal->len);
466
467
468 for(i=0;i<agm->len;i++) {
469 if( should_show[i] == 0 ) {
470 continue;
471 }
472
473 if( i > 0 && should_show[i-1] == 0 ) {
474 /* just had a discontinuity */
475 fprintf(ofp,"\n ---\n");
476 }
477
478 fill_forward_AlignGeneColumn(store,sal,i);
479
480 fprintf(ofp,"%c %4d %2.4f Bits: %2.4f [%1.6f vs %1.6f cum % 8.1f] %2.4f %2.4f | %2.4f %2.4f %2.4f ",
481 should_show[i] == 2 ? '*' : '-',
482 i,
483 agm->forward_coding[i],
484 Score2Bits(Probability2Score(agm->forward_coding[i])),
485 weighted_coding_AlignGeneColumn(store,agmp),
486 weighted_non_coding_AlignGeneColumn(store,agmp),
487 gr == NULL ? 0.0 : Score2Bits(phase_array[i]),
488 agm->splice5_forward[i],
489 agm->splice3_forward[i],
490 agm->reverse_coding[i],agm->splice5_reverse[i],
491 agm->splice3_reverse[i]
492 );
493 for(j=0;j<sal->len;j++) {
494 fprintf(ofp,"%c",sal->seq[j]->seq[i]);
495 }
496 fprintf(ofp,"|");
497
498 for(j=0;j<sal->len;j++) {
499 fprintf(ofp,"%c",char_complement_base(sal->seq[j]->seq[i]));
500 }
501
502 fprintf(ofp," ");
503 if( i > 3 ) {
504 for(j=0;j<sal->len;j++) {
505 if( sal->seq[j]->seq[i-2] == '~' || sal->seq[j]->seq[i-2] == '-' ||
506 sal->seq[j]->seq[i-1] == '~' || sal->seq[j]->seq[i-1] == '-' ||
507 sal->seq[j]->seq[i] == '~' || sal->seq[j]->seq[i] == '-' ) {
508 fprintf(ofp,"?");
509 } else {
510 if( sal->seq[0]->seq[i-2] != sal->seq[1]->seq[i-2] ||
511 sal->seq[0]->seq[i-1] != sal->seq[1]->seq[i-1] ||
512 sal->seq[0]->seq[i] != sal->seq[1]->seq[i]) {
513 fprintf(ofp,"%c",aminoacid_from_seq(ct,sal->seq[j]->seq+i-2));
514 } else {
515 fprintf(ofp,"%c",tolower(aminoacid_from_seq(ct,sal->seq[j]->seq+i-2)));
516 }
517 }
518 }
519 if( agm->forward_coding[i] > 15.0 ) {
520 fprintf(ofp,"*");
521 }
522 }
523 fprintf(ofp,"\n");
524 }
525 }
526
527 # line 582 "aligngenemodel.dy"
create_AlignGeneModel(SeqAlign * sal,AlignGeneModelParam * agmp)528 AlignGeneModel * create_AlignGeneModel(SeqAlign * sal,AlignGeneModelParam * agmp)
529 {
530 AlignGeneModel * out;
531 int i;
532 int j;
533 Probability ss5p;
534 Probability ss3p;
535 AlignGeneColumnStore * store;
536 Sequence * rev;
537
538 assert(agmp);
539 assert(agmp->ct);
540 assert(agmp->protein);
541 assert(agmp->dm);
542
543
544 out = new_AlignGeneModel(sal->seq[0]->len);
545
546 reweight_SeqAlign(sal);
547
548 for(i=0;i<sal->len;i++) {
549 fprintf(stdout,"%s weight %f\n",sal->seq[i]->name,sal->seq[i]->weight);
550 }
551
552 store = new_empty_AlignGeneColumnStore(sal->len);
553
554 for(i=0;i<sal->seq[0]->len;i++) {
555 if( i > 2 ) {
556
557 fill_forward_AlignGeneColumn(store,sal,i);
558
559 /*
560 fprintf(stdout,"%d %f vs %f [%f] | %f vs %f %c%c%c %c%c%c\n",i,
561 simple_coding_AlignGeneColumn(store,agmp),
562 weighted_coding_AlignGeneColumn(store,agmp),
563
564 simple_coding_AlignGeneColumn(store,agmp)/weighted_coding_AlignGeneColumn(store,agmp),
565
566 simple_non_coding_AlignGeneColumn(store,agmp),
567 non_coding_probability_AlignGeneModel(sal,i,agmp),
568 store->codon[0].seq[0],
569 store->codon[0].seq[1],
570 store->codon[0].seq[2],
571 store->codon[1].seq[0],
572 store->codon[1].seq[1],
573 store->codon[1].seq[2]);
574
575
576 */
577
578 out->forward_coding[i] =
579 weighted_coding_AlignGeneColumn(store,agmp)*agmp->total_weight/
580 (weighted_non_coding_AlignGeneColumn(store,agmp));
581
582 }
583
584 ss5p = 1.0;
585 ss3p = 1.0;
586 for(j=0;j<1;j++) {
587 if( sal->seq[j]->seq[i] == 'G' && sal->seq[j]->seq[i+1] == 'T' && i < sal->seq[j]->len-2)
588 ss5p *= prob_SpliceSiteProb(agmp->ss5,sal->seq[j],i);
589 else
590 ss5p = 0.00000001;
591
592 if( sal->seq[j]->seq[i] == 'G' && sal->seq[j]->seq[i-1] == 'A' && i > 2 )
593 ss3p *= prob_SpliceSiteProb(agmp->ss3,sal->seq[j],i);
594 else
595 ss3p = 0.000000001;
596 }
597 if( ss5p > 4 ) {
598 ss5p = 4.0;
599 }
600 if( ss3p > 4 ) {
601 ss3p = 4.0;
602 }
603
604 out->splice5_forward[i] = ss5p;
605 out->splice3_forward[i] = ss3p;
606 }
607
608 rev = reverse_complement_Sequence(sal->seq[0]);
609
610 for(i=0;i<sal->seq[0]->len;i++) {
611 if( i > 2 ) {
612 fill_reverse_AlignGeneColumn(store,sal,i);
613
614 out->reverse_coding[i] =
615 weighted_coding_AlignGeneColumn(store,agmp)*agmp->total_weight/
616 (weighted_non_coding_AlignGeneColumn(store,agmp));
617 }
618
619 if( rev->seq[rev->len-i] == 'G' && rev->seq[rev->len-i-1] == 'T' && i > 2)
620 out->splice5_reverse[i] = prob_SpliceSiteProb(agmp->ss5,rev,rev->len-i);
621 else
622 out->splice5_reverse[i] = 0.00000001;
623
624 if( rev->seq[rev->len-i] == 'G' && rev->seq[rev->len-i+1] == 'A' && i < rev->len -2)
625 out->splice3_reverse[i] = prob_SpliceSiteProb(agmp->ss3,rev,rev->len-i);
626 else
627 out->splice3_reverse[i] = 0.00000001;
628
629
630 if( out->splice5_reverse[i] > 4 ) {
631 out->splice5_reverse[i] = 4.0;
632 }
633
634 if( out->splice3_reverse[i] > 4 ) {
635 out->splice3_reverse[i] = 4.0;
636 }
637
638 /*
639 out->splice5_reverse[i] = 1.0;
640 out->splice3_reverse[i] = 1.0;
641 */
642 }
643
644 free_AlignGeneColumnStore(store);
645 free_Sequence(rev);
646
647 window_AlignGeneModel(sal,out,agmp);
648
649 out->align = hard_link_SeqAlign(sal);
650 out->anchor = hard_link_Sequence(sal->seq[0]);
651
652 return out;
653 }
654
655 # line 709 "aligngenemodel.dy"
prob_SpliceSiteProb(SpliceSiteProb * ssp,Sequence * seq,int pos)656 Probability prob_SpliceSiteProb(SpliceSiteProb * ssp,Sequence * seq,int pos)
657 {
658 assert(ssp);
659 assert(seq);
660
661 if( pos<= ssp->offset ) {
662 return 0.0;
663 }
664
665 return prob_pwmDNA_string(ssp->pwm,seq->seq-ssp->offset+1+pos);
666 }
667
668
669
670
671 # line 724 "aligngenemodel.dy"
non_coding_probability_AlignGeneModel(SeqAlign * sal,int codon_end_pos,AlignGeneModelParam * agmp)672 Probability non_coding_probability_AlignGeneModel(SeqAlign * sal,int codon_end_pos,AlignGeneModelParam * agmp)
673 {
674 int i;
675 int j;
676 int k;
677 Probability out = 1.0;
678 Probability col = 0.0;
679 Probability c = 0.0;
680
681
682 for(k=0;k<3;k++) {
683 for(j=0;j<4;j++) {
684 c= 0.25;
685 for(i=0;i<sal->len;i++) {
686 if( sal->seq[i]->seq[codon_end_pos-k] == '-' ) {
687 return 0.5;
688 }
689 c *= agmp->dm->prob[j][base_from_char(sal->seq[i]->seq[codon_end_pos-k])];
690 }
691 col += c;
692 }
693
694 out *= col;
695 }
696
697 return out;
698
699 }
700
701 # line 753 "aligngenemodel.dy"
coding_probability_AlignGeneModel(SeqAlign * sal,int codon_end_pos,AlignGeneModelParam * agmp)702 Probability coding_probability_AlignGeneModel(SeqAlign * sal,int codon_end_pos,AlignGeneModelParam * agmp)
703 {
704 int i;
705 int j;
706 int codon;
707 Probability out = 0.0;
708 Probability c = 0.0;
709 Probability f;
710
711 assert(sal);
712 assert(agmp);
713 assert(agmp->rm);
714 assert(agmp->protein);
715
716 /*
717 fprintf(stdout,"Testing %c%c%c [%c] vs %c%c%c [%c]\n",
718 sal->seq[0]->seq[codon_end_pos-2],
719 sal->seq[0]->seq[codon_end_pos-1],
720 sal->seq[0]->seq[codon_end_pos],
721 aminoacid_from_codon(agmp->ct,codon_from_seq(sal->seq[0]->seq+codon_end_pos-2)),
722 sal->seq[1]->seq[codon_end_pos-2],
723 sal->seq[1]->seq[codon_end_pos-1],
724 sal->seq[1]->seq[codon_end_pos],
725 aminoacid_from_codon(agmp->ct,codon_from_seq(sal->seq[1]->seq+codon_end_pos-2)) );
726 */
727
728 for(j=0;j<26;j++) {
729 if( agmp->rm->aminoacid[j] < 0.0000000001 ) {
730 continue;
731 }
732 c = agmp->rm->aminoacid[j];
733 for(i=0;i<sal->len;i++) {
734 if( !isalpha(sal->seq[i]->seq[codon_end_pos-2]) ||
735 !isalpha(sal->seq[i]->seq[codon_end_pos-1]) ||
736 !isalpha(sal->seq[i]->seq[codon_end_pos]) ) {
737 return 0.0;
738 }
739
740 codon = codon_from_seq(sal->seq[i]->seq+codon_end_pos-2);
741 if( is_stop_codon(codon,agmp->ct) ) {
742 return 0.0;
743 }
744 f = agmp->protein->comp[j][aminoacid_from_codon(agmp->ct,codon)-'A'];
745
746 /*
747 fprintf(stdout,"For amino acid %c, match %.4f cum %.4f\n",j+'A',
748 agmp->protein->comp[j][aminoacid_from_codon(agmp->ct,codon)-'A'],
749 c);
750 */
751
752 c = c*f;
753 }
754 out += c;
755 }
756
757 /* fprintf(stdout,"final %.4f\n",out);*/
758 return out;
759
760 }
761
762
763 # line 814 "aligngenemodel.dy"
new_AlignGeneModelParam_from_argv(int * argc,char ** argv)764 AlignGeneModelParam * new_AlignGeneModelParam_from_argv(int * argc,char ** argv)
765 {
766 AlignGeneModelParam * out;
767 CompProb * cp;
768 DnaProbMatrix * dm;
769
770 GeneModelParam * gmp;
771 GeneStats * gs;
772
773 char * temp;
774 double temp_prob;
775
776 out = AlignGeneModelParam_alloc();
777
778 if( (temp = strip_out_assigned_argument(argc,argv,"am_codon")) == NULL ) {
779 temp = "codon.table";
780 }
781 out->ct = read_CodonTable_file(temp);
782 assert(out->ct);
783
784 if( (temp = strip_out_assigned_argument(argc,argv,"am_protein")) == NULL ) {
785 temp = "wag85";
786 }
787 out->protein = read_Blast_file_CompProb(temp);
788 assert(out->protein);
789
790 temp_prob = 0.85;
791 strip_out_float_argument(argc,argv,"am_dna",&temp_prob);
792
793 out->dm = DnaProbMatrix_from_match(0.85,NMaskType_VARIABLE);
794
795 gmp = new_GeneModelParam_from_argv(argc,argv);
796
797 if((gs=GeneStats_from_GeneModelParam(gmp)) == NULL ) {
798 fatal("Could not build gene stats");
799 }
800
801 out->gs = gs;
802 out->rm = default_RandomModel();
803
804 out->ss5 = SpliceSiteProb_alloc();
805 out->ss5->pwm = pwmDNA_from_SeqAlign(gs->splice5,1);
806
807 out->ss5->offset = gs->splice5_offset;
808 fold_randommodel_pwmDNA(out->ss5->pwm,gs->rnd);
809
810 out->ss3 = SpliceSiteProb_alloc();
811 out->ss3->pwm = pwmDNA_from_SeqAlign(gs->splice3,1);
812
813 out->ss3->offset = gs->splice3_offset;
814 fold_randommodel_pwmDNA(out->ss3->pwm,gs->rnd);
815
816
817 out->nonanchor_stop = 0.001;
818 strip_out_float_argument(argc,argv,"am_nonanchor_stop_pen",&out->nonanchor_stop);
819
820 out->tolerate_nonanchor_stops = 1;
821 strip_out_boolean_def_argument(argc,argv,"am_nonanchor_stop",&out->tolerate_nonanchor_stops);
822
823 out->coding_window_thres = 3.0;
824 strip_out_float_argument(argc,argv,"am_coding_bonus_thres",&out->coding_window_thres);
825
826
827 out->noncoding_window_thres = 2.0;
828 strip_out_float_argument(argc,argv,"am_noncoding_bonus_thres",&out->noncoding_window_thres);
829
830
831 out->coding_window_bonus = 4.0;
832 strip_out_float_argument(argc,argv,"am_coding_bonus_factor",&out->coding_window_bonus);
833
834 out->noncoding_window_pen = 0.25;
835 strip_out_float_argument(argc,argv,"am_noncoding_bonus_factor",&out->noncoding_window_pen);
836
837
838 out->total_weight = 0.5;
839 strip_out_float_argument(argc,argv,"am_overall_coding",&out->total_weight);
840
841
842 return out;
843
844 }
845
846 # line 896 "aligngenemodel.dy"
show_help_AlignGeneModelParam(FILE * ofp)847 void show_help_AlignGeneModelParam(FILE * ofp)
848 {
849 fprintf(ofp,"Align Gene Model Parameters\n");
850 fprintf(ofp," -am_codon [codon.table] codon table\n");
851 fprintf(ofp," -am_protein [wag85] protein comparison probabilities\n");
852 fprintf(ofp," -am_dna [0.85] DNA non coding model match probability\n");
853 fprintf(ofp," -[no]am_nonanchor_stop allow stops in non anchor sequences [true]\n");
854 fprintf(ofp," -am_nonanchor_stop_pen penalty for non anchor stops when allowed [0.001]\n");
855 fprintf(ofp," -am_coding_bonus_thres coding threshold to use conservation window [3.0]\n");
856 fprintf(ofp," -am_noncoding_bonus_thres nonc threshold to use conservation window [2.0]\n");
857 fprintf(ofp," -am_coding_bonus_factor coding bonus factor when windowed [4.0]\n");
858 fprintf(ofp," -am_noncoding_bonus_factor nonc bonus factor when windowed [0.25]\n");
859 fprintf(ofp," -am_overall_coding overall per column \"odds-prior\" for coding [0.5]\n");
860
861 fprintf(ofp,"Gene Model statistics used by AlignModel\n");
862 show_help_GeneModelParam(ofp);
863
864 return;
865
866 }
867
868
869 # line 918 "aligngenemodel.dy"
std_AlignGeneModelParam(CompProb * cp,DnaProbMatrix * dm,CodonTable * ct,GeneStats * gs,Probability change)870 AlignGeneModelParam * std_AlignGeneModelParam(CompProb * cp,DnaProbMatrix * dm,CodonTable * ct,GeneStats * gs,Probability change)
871 {
872 AlignGeneModelParam * out;
873
874 out = AlignGeneModelParam_alloc();
875
876 assert(cp);
877 assert(dm);
878 assert(ct);
879 assert(gs);
880
881 assert(gs->splice5);
882 assert(gs->splice5->len > 1);
883 assert(gs->splice5->seq[0]->len > 0);
884
885 assert(gs->splice3);
886 assert(gs->splice3->len > 1);
887 assert(gs->splice3->seq[0]->len > 0);
888
889 out->ncsm = create_NonCodingSimpleModel(change);
890 out->protein = hard_link_CompProb(cp);
891 out->rm = default_RandomModel();
892 out->ct = hard_link_CodonTable(ct);
893 out->dm = hard_link_DnaProbMatrix(dm);
894
895 out->ss5 = SpliceSiteProb_alloc();
896 out->ss5->pwm = pwmDNA_from_SeqAlign(gs->splice5,1);
897
898 out->ss5->offset = gs->splice5_offset;
899 fold_randommodel_pwmDNA(out->ss5->pwm,gs->rnd);
900
901 /* show_pwmDNA_col(out->ss5->pwm,stdout);
902 fprintf(stdout,"\n");
903 */
904 out->ss3 = SpliceSiteProb_alloc();
905 out->ss3->pwm = pwmDNA_from_SeqAlign(gs->splice3,1);
906
907
908
909 out->ss3->offset = gs->splice3_offset;
910 fold_randommodel_pwmDNA(out->ss3->pwm,gs->rnd);
911 /*
912 show_pwmDNA_col(out->ss3->pwm,stdout);
913 fprintf(stdout,"\n");
914 */
915
916 out->nonanchor_stop = 0.001;
917 out->coding_window_thres = 3.0;
918 out->noncoding_window_thres = 2.0;
919 out->coding_window_bonus = 4.0;
920 out->noncoding_window_pen = 0.25;
921 out->tolerate_nonanchor_stops = 1;
922 out->total_weight = 0.5;
923
924 return out;
925 }
926
927 # line 975 "aligngenemodel.dy"
AlignGeneModelScore_from_AlignGeneModel(AlignGeneModel * agm)928 AlignGeneModelScore * AlignGeneModelScore_from_AlignGeneModel(AlignGeneModel * agm)
929 {
930 AlignGeneModelScore * out;
931
932 out = AlignGeneModelScore_alloc();
933
934 out->len = agm->len;
935
936 out->forward_coding = calloc(agm->len,sizeof(Score));
937 out->reverse_coding = calloc(agm->len,sizeof(Score));
938 out->splice5_forward = calloc(agm->len,sizeof(Score));
939 out->splice3_forward = calloc(agm->len,sizeof(Score));
940 out->splice5_reverse = calloc(agm->len,sizeof(Score));
941 out->splice3_reverse = calloc(agm->len,sizeof(Score));
942
943 Probability2Score_move(agm->forward_coding,out->forward_coding,agm->len);
944 Probability2Score_move(agm->reverse_coding,out->reverse_coding,agm->len);
945 Probability2Score_move(agm->splice5_forward,out->splice5_forward,agm->len);
946 Probability2Score_move(agm->splice3_forward,out->splice3_forward,agm->len);
947 Probability2Score_move(agm->splice5_reverse,out->splice5_reverse,agm->len);
948 Probability2Score_move(agm->splice3_reverse,out->splice3_reverse,agm->len);
949
950
951 out->align = hard_link_SeqAlign(agm->align);
952 out->anchor = hard_link_Sequence(agm->anchor);
953
954 return out;
955 }
956
957 # line 1004 "aligngenemodel.dy"
free_Probability(Probability * p)958 Probability * free_Probability(Probability * p)
959 {
960 ckfree(p);
961 return NULL;
962 }
963
964 # line 1010 "aligngenemodel.dy"
free_Score(Score * p)965 Score * free_Score(Score * p)
966 {
967 ckfree(p);
968 return NULL;
969 }
970 # line 969 "aligngenemodel.c"
971 /* Function: hard_link_NonCodingSimpleModel(obj)
972 *
973 * Descrip: Bumps up the reference count of the object
974 * Meaning that multiple pointers can 'own' it
975 *
976 *
977 * Arg: obj [UNKN ] Object to be hard linked [NonCodingSimpleModel *]
978 *
979 * Return [UNKN ] Undocumented return value [NonCodingSimpleModel *]
980 *
981 */
hard_link_NonCodingSimpleModel(NonCodingSimpleModel * obj)982 NonCodingSimpleModel * hard_link_NonCodingSimpleModel(NonCodingSimpleModel * obj)
983 {
984 if( obj == NULL ) {
985 warn("Trying to hard link to a NonCodingSimpleModel object: passed a NULL object");
986 return NULL;
987 }
988 obj->dynamite_hard_link++;
989 return obj;
990 }
991
992
993 /* Function: NonCodingSimpleModel_alloc(void)
994 *
995 * Descrip: Allocates structure: assigns defaults if given
996 *
997 *
998 *
999 * Return [UNKN ] Undocumented return value [NonCodingSimpleModel *]
1000 *
1001 */
NonCodingSimpleModel_alloc(void)1002 NonCodingSimpleModel * NonCodingSimpleModel_alloc(void)
1003 {
1004 NonCodingSimpleModel * out; /* out is exported at end of function */
1005
1006
1007 /* call ckalloc and see if NULL */
1008 if((out=(NonCodingSimpleModel *) ckalloc (sizeof(NonCodingSimpleModel))) == NULL) {
1009 warn("NonCodingSimpleModel_alloc failed ");
1010 return NULL; /* calling function should respond! */
1011 }
1012 out->dynamite_hard_link = 1;
1013 #ifdef PTHREAD
1014 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1015 #endif
1016 out->identical = 0.0;
1017 out->one_off = 0.0;
1018 out->two_off = 0.0;
1019
1020
1021 return out;
1022 }
1023
1024
1025 /* Function: free_NonCodingSimpleModel(obj)
1026 *
1027 * Descrip: Free Function: removes the memory held by obj
1028 * Will chain up to owned members and clear all lists
1029 *
1030 *
1031 * Arg: obj [UNKN ] Object that is free'd [NonCodingSimpleModel *]
1032 *
1033 * Return [UNKN ] Undocumented return value [NonCodingSimpleModel *]
1034 *
1035 */
free_NonCodingSimpleModel(NonCodingSimpleModel * obj)1036 NonCodingSimpleModel * free_NonCodingSimpleModel(NonCodingSimpleModel * obj)
1037 {
1038 int return_early = 0;
1039
1040
1041 if( obj == NULL) {
1042 warn("Attempting to free a NULL pointer to a NonCodingSimpleModel obj. Should be trappable");
1043 return NULL;
1044 }
1045
1046
1047 #ifdef PTHREAD
1048 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1049 #endif
1050 if( obj->dynamite_hard_link > 1) {
1051 return_early = 1;
1052 obj->dynamite_hard_link--;
1053 }
1054 #ifdef PTHREAD
1055 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1056 #endif
1057 if( return_early == 1)
1058 return NULL;
1059
1060
1061 ckfree(obj);
1062 return NULL;
1063 }
1064
1065
1066 /* Function: hard_link_SpliceSiteProb(obj)
1067 *
1068 * Descrip: Bumps up the reference count of the object
1069 * Meaning that multiple pointers can 'own' it
1070 *
1071 *
1072 * Arg: obj [UNKN ] Object to be hard linked [SpliceSiteProb *]
1073 *
1074 * Return [UNKN ] Undocumented return value [SpliceSiteProb *]
1075 *
1076 */
hard_link_SpliceSiteProb(SpliceSiteProb * obj)1077 SpliceSiteProb * hard_link_SpliceSiteProb(SpliceSiteProb * obj)
1078 {
1079 if( obj == NULL ) {
1080 warn("Trying to hard link to a SpliceSiteProb object: passed a NULL object");
1081 return NULL;
1082 }
1083 obj->dynamite_hard_link++;
1084 return obj;
1085 }
1086
1087
1088 /* Function: SpliceSiteProb_alloc(void)
1089 *
1090 * Descrip: Allocates structure: assigns defaults if given
1091 *
1092 *
1093 *
1094 * Return [UNKN ] Undocumented return value [SpliceSiteProb *]
1095 *
1096 */
SpliceSiteProb_alloc(void)1097 SpliceSiteProb * SpliceSiteProb_alloc(void)
1098 {
1099 SpliceSiteProb * out; /* out is exported at end of function */
1100
1101
1102 /* call ckalloc and see if NULL */
1103 if((out=(SpliceSiteProb *) ckalloc (sizeof(SpliceSiteProb))) == NULL) {
1104 warn("SpliceSiteProb_alloc failed ");
1105 return NULL; /* calling function should respond! */
1106 }
1107 out->dynamite_hard_link = 1;
1108 #ifdef PTHREAD
1109 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1110 #endif
1111 out->pwm = NULL;
1112 out->offset = 0;
1113
1114
1115 return out;
1116 }
1117
1118
1119 /* Function: free_SpliceSiteProb(obj)
1120 *
1121 * Descrip: Free Function: removes the memory held by obj
1122 * Will chain up to owned members and clear all lists
1123 *
1124 *
1125 * Arg: obj [UNKN ] Object that is free'd [SpliceSiteProb *]
1126 *
1127 * Return [UNKN ] Undocumented return value [SpliceSiteProb *]
1128 *
1129 */
free_SpliceSiteProb(SpliceSiteProb * obj)1130 SpliceSiteProb * free_SpliceSiteProb(SpliceSiteProb * obj)
1131 {
1132 int return_early = 0;
1133
1134
1135 if( obj == NULL) {
1136 warn("Attempting to free a NULL pointer to a SpliceSiteProb obj. Should be trappable");
1137 return NULL;
1138 }
1139
1140
1141 #ifdef PTHREAD
1142 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1143 #endif
1144 if( obj->dynamite_hard_link > 1) {
1145 return_early = 1;
1146 obj->dynamite_hard_link--;
1147 }
1148 #ifdef PTHREAD
1149 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1150 #endif
1151 if( return_early == 1)
1152 return NULL;
1153 if( obj->pwm != NULL)
1154 free_pwmDNA(obj->pwm);
1155
1156
1157 ckfree(obj);
1158 return NULL;
1159 }
1160
1161
1162 /* Function: hard_link_AlignGeneColumnStore(obj)
1163 *
1164 * Descrip: Bumps up the reference count of the object
1165 * Meaning that multiple pointers can 'own' it
1166 *
1167 *
1168 * Arg: obj [UNKN ] Object to be hard linked [AlignGeneColumnStore *]
1169 *
1170 * Return [UNKN ] Undocumented return value [AlignGeneColumnStore *]
1171 *
1172 */
hard_link_AlignGeneColumnStore(AlignGeneColumnStore * obj)1173 AlignGeneColumnStore * hard_link_AlignGeneColumnStore(AlignGeneColumnStore * obj)
1174 {
1175 if( obj == NULL ) {
1176 warn("Trying to hard link to a AlignGeneColumnStore object: passed a NULL object");
1177 return NULL;
1178 }
1179 obj->dynamite_hard_link++;
1180 return obj;
1181 }
1182
1183
1184 /* Function: AlignGeneColumnStore_alloc(void)
1185 *
1186 * Descrip: Allocates structure: assigns defaults if given
1187 *
1188 *
1189 *
1190 * Return [UNKN ] Undocumented return value [AlignGeneColumnStore *]
1191 *
1192 */
AlignGeneColumnStore_alloc(void)1193 AlignGeneColumnStore * AlignGeneColumnStore_alloc(void)
1194 {
1195 AlignGeneColumnStore * out; /* out is exported at end of function */
1196
1197
1198 /* call ckalloc and see if NULL */
1199 if((out=(AlignGeneColumnStore *) ckalloc (sizeof(AlignGeneColumnStore))) == NULL) {
1200 warn("AlignGeneColumnStore_alloc failed ");
1201 return NULL; /* calling function should respond! */
1202 }
1203 out->dynamite_hard_link = 1;
1204 #ifdef PTHREAD
1205 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1206 #endif
1207 out->codon = NULL;
1208 out->len = 0;
1209
1210
1211 return out;
1212 }
1213
1214
1215 /* Function: free_AlignGeneColumnStore(obj)
1216 *
1217 * Descrip: Free Function: removes the memory held by obj
1218 * Will chain up to owned members and clear all lists
1219 *
1220 *
1221 * Arg: obj [UNKN ] Object that is free'd [AlignGeneColumnStore *]
1222 *
1223 * Return [UNKN ] Undocumented return value [AlignGeneColumnStore *]
1224 *
1225 */
free_AlignGeneColumnStore(AlignGeneColumnStore * obj)1226 AlignGeneColumnStore * free_AlignGeneColumnStore(AlignGeneColumnStore * obj)
1227 {
1228 int return_early = 0;
1229
1230
1231 if( obj == NULL) {
1232 warn("Attempting to free a NULL pointer to a AlignGeneColumnStore obj. Should be trappable");
1233 return NULL;
1234 }
1235
1236
1237 #ifdef PTHREAD
1238 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1239 #endif
1240 if( obj->dynamite_hard_link > 1) {
1241 return_early = 1;
1242 obj->dynamite_hard_link--;
1243 }
1244 #ifdef PTHREAD
1245 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1246 #endif
1247 if( return_early == 1)
1248 return NULL;
1249 if( obj->codon != NULL)
1250 free_AlignGeneCodon(obj->codon);
1251
1252
1253 ckfree(obj);
1254 return NULL;
1255 }
1256
1257
1258 /* Function: hard_link_AlignGeneModel(obj)
1259 *
1260 * Descrip: Bumps up the reference count of the object
1261 * Meaning that multiple pointers can 'own' it
1262 *
1263 *
1264 * Arg: obj [UNKN ] Object to be hard linked [AlignGeneModel *]
1265 *
1266 * Return [UNKN ] Undocumented return value [AlignGeneModel *]
1267 *
1268 */
hard_link_AlignGeneModel(AlignGeneModel * obj)1269 AlignGeneModel * hard_link_AlignGeneModel(AlignGeneModel * obj)
1270 {
1271 if( obj == NULL ) {
1272 warn("Trying to hard link to a AlignGeneModel object: passed a NULL object");
1273 return NULL;
1274 }
1275 obj->dynamite_hard_link++;
1276 return obj;
1277 }
1278
1279
1280 /* Function: AlignGeneModel_alloc(void)
1281 *
1282 * Descrip: Allocates structure: assigns defaults if given
1283 *
1284 *
1285 *
1286 * Return [UNKN ] Undocumented return value [AlignGeneModel *]
1287 *
1288 */
AlignGeneModel_alloc(void)1289 AlignGeneModel * AlignGeneModel_alloc(void)
1290 {
1291 AlignGeneModel * out; /* out is exported at end of function */
1292
1293
1294 /* call ckalloc and see if NULL */
1295 if((out=(AlignGeneModel *) ckalloc (sizeof(AlignGeneModel))) == NULL) {
1296 warn("AlignGeneModel_alloc failed ");
1297 return NULL; /* calling function should respond! */
1298 }
1299 out->dynamite_hard_link = 1;
1300 #ifdef PTHREAD
1301 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1302 #endif
1303 out->len = 0;
1304 out->forward_coding = NULL;
1305 out->reverse_coding = NULL;
1306 out->splice5_forward = NULL;
1307 out->splice3_forward = NULL;
1308 out->splice5_reverse = NULL;
1309 out->splice3_reverse = NULL;
1310 out->align = NULL;
1311 out->anchor = NULL;
1312 out->change_rate = NULL;
1313
1314
1315 return out;
1316 }
1317
1318
1319 /* Function: free_AlignGeneModel(obj)
1320 *
1321 * Descrip: Free Function: removes the memory held by obj
1322 * Will chain up to owned members and clear all lists
1323 *
1324 *
1325 * Arg: obj [UNKN ] Object that is free'd [AlignGeneModel *]
1326 *
1327 * Return [UNKN ] Undocumented return value [AlignGeneModel *]
1328 *
1329 */
free_AlignGeneModel(AlignGeneModel * obj)1330 AlignGeneModel * free_AlignGeneModel(AlignGeneModel * obj)
1331 {
1332 int return_early = 0;
1333
1334
1335 if( obj == NULL) {
1336 warn("Attempting to free a NULL pointer to a AlignGeneModel obj. Should be trappable");
1337 return NULL;
1338 }
1339
1340
1341 #ifdef PTHREAD
1342 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1343 #endif
1344 if( obj->dynamite_hard_link > 1) {
1345 return_early = 1;
1346 obj->dynamite_hard_link--;
1347 }
1348 #ifdef PTHREAD
1349 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1350 #endif
1351 if( return_early == 1)
1352 return NULL;
1353 if( obj->forward_coding != NULL)
1354 free_Probability(obj->forward_coding);
1355 if( obj->reverse_coding != NULL)
1356 free_Probability(obj->reverse_coding);
1357 if( obj->splice5_forward != NULL)
1358 free_Probability(obj->splice5_forward);
1359 if( obj->splice3_forward != NULL)
1360 free_Probability(obj->splice3_forward);
1361 if( obj->splice5_reverse != NULL)
1362 free_Probability(obj->splice5_reverse);
1363 if( obj->splice3_reverse != NULL)
1364 free_Probability(obj->splice3_reverse);
1365 if( obj->align != NULL)
1366 free_SeqAlign(obj->align);
1367 if( obj->anchor != NULL)
1368 free_Sequence(obj->anchor);
1369 if( obj->change_rate != NULL)
1370 ckfree(obj->change_rate);
1371
1372
1373 ckfree(obj);
1374 return NULL;
1375 }
1376
1377
1378 /* Function: hard_link_AlignGeneModelScore(obj)
1379 *
1380 * Descrip: Bumps up the reference count of the object
1381 * Meaning that multiple pointers can 'own' it
1382 *
1383 *
1384 * Arg: obj [UNKN ] Object to be hard linked [AlignGeneModelScore *]
1385 *
1386 * Return [UNKN ] Undocumented return value [AlignGeneModelScore *]
1387 *
1388 */
hard_link_AlignGeneModelScore(AlignGeneModelScore * obj)1389 AlignGeneModelScore * hard_link_AlignGeneModelScore(AlignGeneModelScore * obj)
1390 {
1391 if( obj == NULL ) {
1392 warn("Trying to hard link to a AlignGeneModelScore object: passed a NULL object");
1393 return NULL;
1394 }
1395 obj->dynamite_hard_link++;
1396 return obj;
1397 }
1398
1399
1400 /* Function: AlignGeneModelScore_alloc(void)
1401 *
1402 * Descrip: Allocates structure: assigns defaults if given
1403 *
1404 *
1405 *
1406 * Return [UNKN ] Undocumented return value [AlignGeneModelScore *]
1407 *
1408 */
AlignGeneModelScore_alloc(void)1409 AlignGeneModelScore * AlignGeneModelScore_alloc(void)
1410 {
1411 AlignGeneModelScore * out; /* out is exported at end of function */
1412
1413
1414 /* call ckalloc and see if NULL */
1415 if((out=(AlignGeneModelScore *) ckalloc (sizeof(AlignGeneModelScore))) == NULL) {
1416 warn("AlignGeneModelScore_alloc failed ");
1417 return NULL; /* calling function should respond! */
1418 }
1419 out->dynamite_hard_link = 1;
1420 #ifdef PTHREAD
1421 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1422 #endif
1423 out->len = 0;
1424 out->forward_coding = NULL;
1425 out->reverse_coding = NULL;
1426 out->splice5_forward = NULL;
1427 out->splice3_forward = NULL;
1428 out->splice5_reverse = NULL;
1429 out->splice3_reverse = NULL;
1430 out->align = NULL;
1431 out->anchor = NULL;
1432
1433
1434 return out;
1435 }
1436
1437
1438 /* Function: free_AlignGeneModelScore(obj)
1439 *
1440 * Descrip: Free Function: removes the memory held by obj
1441 * Will chain up to owned members and clear all lists
1442 *
1443 *
1444 * Arg: obj [UNKN ] Object that is free'd [AlignGeneModelScore *]
1445 *
1446 * Return [UNKN ] Undocumented return value [AlignGeneModelScore *]
1447 *
1448 */
free_AlignGeneModelScore(AlignGeneModelScore * obj)1449 AlignGeneModelScore * free_AlignGeneModelScore(AlignGeneModelScore * obj)
1450 {
1451 int return_early = 0;
1452
1453
1454 if( obj == NULL) {
1455 warn("Attempting to free a NULL pointer to a AlignGeneModelScore obj. Should be trappable");
1456 return NULL;
1457 }
1458
1459
1460 #ifdef PTHREAD
1461 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1462 #endif
1463 if( obj->dynamite_hard_link > 1) {
1464 return_early = 1;
1465 obj->dynamite_hard_link--;
1466 }
1467 #ifdef PTHREAD
1468 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1469 #endif
1470 if( return_early == 1)
1471 return NULL;
1472 if( obj->forward_coding != NULL)
1473 free_Score(obj->forward_coding);
1474 if( obj->reverse_coding != NULL)
1475 free_Score(obj->reverse_coding);
1476 if( obj->splice5_forward != NULL)
1477 free_Score(obj->splice5_forward);
1478 if( obj->splice3_forward != NULL)
1479 free_Score(obj->splice3_forward);
1480 if( obj->splice5_reverse != NULL)
1481 free_Score(obj->splice5_reverse);
1482 if( obj->splice3_reverse != NULL)
1483 free_Score(obj->splice3_reverse);
1484 if( obj->align != NULL)
1485 free_SeqAlign(obj->align);
1486 if( obj->anchor != NULL)
1487 free_Sequence(obj->anchor);
1488
1489
1490 ckfree(obj);
1491 return NULL;
1492 }
1493
1494
1495 /* Function: hard_link_AlignGeneModelParam(obj)
1496 *
1497 * Descrip: Bumps up the reference count of the object
1498 * Meaning that multiple pointers can 'own' it
1499 *
1500 *
1501 * Arg: obj [UNKN ] Object to be hard linked [AlignGeneModelParam *]
1502 *
1503 * Return [UNKN ] Undocumented return value [AlignGeneModelParam *]
1504 *
1505 */
hard_link_AlignGeneModelParam(AlignGeneModelParam * obj)1506 AlignGeneModelParam * hard_link_AlignGeneModelParam(AlignGeneModelParam * obj)
1507 {
1508 if( obj == NULL ) {
1509 warn("Trying to hard link to a AlignGeneModelParam object: passed a NULL object");
1510 return NULL;
1511 }
1512 obj->dynamite_hard_link++;
1513 return obj;
1514 }
1515
1516
1517 /* Function: AlignGeneModelParam_alloc(void)
1518 *
1519 * Descrip: Allocates structure: assigns defaults if given
1520 *
1521 *
1522 *
1523 * Return [UNKN ] Undocumented return value [AlignGeneModelParam *]
1524 *
1525 */
AlignGeneModelParam_alloc(void)1526 AlignGeneModelParam * AlignGeneModelParam_alloc(void)
1527 {
1528 AlignGeneModelParam * out; /* out is exported at end of function */
1529
1530
1531 /* call ckalloc and see if NULL */
1532 if((out=(AlignGeneModelParam *) ckalloc (sizeof(AlignGeneModelParam))) == NULL) {
1533 warn("AlignGeneModelParam_alloc failed ");
1534 return NULL; /* calling function should respond! */
1535 }
1536 out->dynamite_hard_link = 1;
1537 #ifdef PTHREAD
1538 pthread_mutex_init(&(out->dynamite_mutex),NULL);
1539 #endif
1540 out->protein = NULL;
1541 out->rm = NULL;
1542 out->cm = NULL;
1543 out->ct = NULL;
1544 out->dm = NULL;
1545 out->ss5 = NULL;
1546 out->ss3 = NULL;
1547 out->gs = NULL;
1548 out->total_weight = 2;
1549 out->tolerate_nonanchor_stops = TRUE;
1550 out->nonanchor_stop = 0.0;
1551 out->ncsm = NULL;
1552 out->coding_window_thres = 0;
1553 out->noncoding_window_thres = 0;
1554 out->coding_window_bonus = 0;
1555 out->noncoding_window_pen = 0;
1556
1557
1558 return out;
1559 }
1560
1561
1562 /* Function: free_AlignGeneModelParam(obj)
1563 *
1564 * Descrip: Free Function: removes the memory held by obj
1565 * Will chain up to owned members and clear all lists
1566 *
1567 *
1568 * Arg: obj [UNKN ] Object that is free'd [AlignGeneModelParam *]
1569 *
1570 * Return [UNKN ] Undocumented return value [AlignGeneModelParam *]
1571 *
1572 */
free_AlignGeneModelParam(AlignGeneModelParam * obj)1573 AlignGeneModelParam * free_AlignGeneModelParam(AlignGeneModelParam * obj)
1574 {
1575 int return_early = 0;
1576
1577
1578 if( obj == NULL) {
1579 warn("Attempting to free a NULL pointer to a AlignGeneModelParam obj. Should be trappable");
1580 return NULL;
1581 }
1582
1583
1584 #ifdef PTHREAD
1585 assert(pthread_mutex_lock(&(obj->dynamite_mutex)) == 0);
1586 #endif
1587 if( obj->dynamite_hard_link > 1) {
1588 return_early = 1;
1589 obj->dynamite_hard_link--;
1590 }
1591 #ifdef PTHREAD
1592 assert(pthread_mutex_unlock(&(obj->dynamite_mutex)) == 0);
1593 #endif
1594 if( return_early == 1)
1595 return NULL;
1596 if( obj->protein != NULL)
1597 free_CompProb(obj->protein);
1598 if( obj->rm != NULL)
1599 free_RandomModel(obj->rm);
1600 if( obj->cm != NULL)
1601 free_CodonMapper(obj->cm);
1602 if( obj->ct != NULL)
1603 free_CodonTable(obj->ct);
1604 if( obj->dm != NULL)
1605 free_DnaProbMatrix(obj->dm);
1606 if( obj->ss5 != NULL)
1607 free_SpliceSiteProb(obj->ss5);
1608 if( obj->ss3 != NULL)
1609 free_SpliceSiteProb(obj->ss3);
1610 if( obj->gs != NULL)
1611 free_GeneStats(obj->gs);
1612 if( obj->ncsm != NULL)
1613 free_NonCodingSimpleModel(obj->ncsm);
1614
1615
1616 ckfree(obj);
1617 return NULL;
1618 }
1619
1620
1621
1622 #ifdef _cplusplus
1623 }
1624 #endif
1625