1 /* @source cpgplot application
2 **
3 ** Plots CpG island areas
4 **
5 ** @author Copyright (C) Alan Bleasby (ableasby@hgmp.mrc.ac.uk)
6 ** @@
7 **
8 ** This program is free software; you can redistribute it and/or
9 ** modify it under the terms of the GNU General Public License
10 ** as published by the Free Software Foundation; either version 2
11 ** of the License, or (at your option) any later version.
12 **
13 ** This program is distributed in the hope that it will be useful,
14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 ** GNU General Public License for more details.
17 **
18 ** You should have received a copy of the GNU General Public License
19 ** along with this program; if not, write to the Free Software
20 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21 ******************************************************************************/
22
23 #include "emboss.h"
24 #include <math.h>
25 #include <stdlib.h>
26
27
28
29
30 static void cpgplot_findbases(const AjPStr substr, ajint len,
31 ajint window, ajint shift, float *obsexp,
32 float *xypc, const AjPStr bases,
33 float *obsexpmax,
34 ajint *plotstart, ajint *plotend);
35
36 static void cpgplot_countbases(const char *seq, const char *bases,
37 ajint window,
38 float *cx, float *cy, float *cxpy);
39
40 static void cpgplot_identify(AjPFile outf,
41 const float *obsexp, const float *xypc,
42 AjBool *thresh, ajint begin, ajint len,
43 ajint shift, const char *bases, const char *name,
44 ajint minlen, float minobsexp, float minpc,
45 AjPFeattabOut featout);
46
47 static void cpgplot_reportislands(AjPFile outf, const AjBool *thresh,
48 const char *bases, const char *name,
49 float minobsexp, float minpc,
50 ajint minlen, ajint begin, ajint len);
51
52 static void cpgplot_plotit(const char *seq,
53 ajint begin, ajint len,
54 const float *obsexp, const float *xypc,
55 const AjBool *thresh,
56 float obsexpmax, AjBool doobsexp,
57 AjBool docg, AjBool dopc, AjPGraph mult);
58
59 static void cpgplot_dumpfeatout(AjPFeattabOut featout, const AjBool *thresh,
60 const char *seqname, ajint begin, ajint len);
61
62
63
64
65 /* @prog cpgplot **************************************************************
66 **
67 ** Plot CpG rich areas
68 **
69 ******************************************************************************/
70
main(int argc,char ** argv)71 int main(int argc, char **argv)
72 {
73 AjPSeqall seqall;
74 AjPSeq seq = NULL;
75 AjPFile outf = NULL;
76 AjPStr strand = NULL;
77 AjPStr substr = NULL;
78 AjPStr bases = NULL;
79 AjPGraph mult;
80 AjBool doobsexp;
81 AjBool dopc;
82 AjBool docg;
83
84 ajint begin;
85 ajint end;
86 ajint len;
87
88 ajint minlen;
89 float minobsexp;
90 float minpc;
91 AjBool doplot;
92
93 ajint window;
94 ajint shift;
95 ajint plotstart;
96 ajint plotend;
97
98 float *xypc = NULL;
99 float *obsexp = NULL;
100 AjBool *thresh = NULL;
101 float obsexpmax;
102
103 ajint i;
104 ajint maxarr;
105 AjPFeattabOut featout = NULL;
106
107
108 embInit("cpgplot",argc,argv);
109
110 seqall = ajAcdGetSeqall("sequence");
111 window = ajAcdGetInt("window");
112 shift = 1; /* other values broken - needs rewrite to fix*/
113 outf = ajAcdGetOutfile("outfile");
114 minobsexp = ajAcdGetFloat("minoe");
115 minlen = ajAcdGetInt("minlen");
116 minpc = ajAcdGetFloat("minpc");
117 doplot = ajAcdGetToggle("plot");
118 mult = ajAcdGetGraphxy ("graph");
119 doobsexp = ajAcdGetBoolean("obsexp");
120 docg = ajAcdGetBoolean("cg");
121 dopc = ajAcdGetBoolean("pc");
122 featout = ajAcdGetFeatout("outfeat");
123
124
125
126 substr = ajStrNew();
127 bases = ajStrNewC("CG");
128
129
130 maxarr = 0;
131
132 while(ajSeqallNext(seqall, &seq))
133 {
134 begin = ajSeqallGetseqBegin(seqall);
135 end = ajSeqallGetseqEnd(seqall);
136 strand = ajSeqGetSeqCopyS(seq);
137 ajStrFmtUpper(&strand);
138
139 ajStrAssignSubC(&substr,ajStrGetPtr(strand),--begin,--end);
140 len = ajStrGetLen(substr);
141
142 if(len > maxarr)
143 {
144 AJCRESIZE(obsexp,len);
145 AJCRESIZE(thresh,len);
146 AJCRESIZE(xypc,len);
147 }
148
149 for(i=0;i<len;++i)
150 {
151 thresh[i] = ajFalse;
152 obsexp[i] = xypc[i] = 0.0;
153 }
154
155
156
157 cpgplot_findbases(substr, len, window, shift, obsexp, xypc,
158 bases, &obsexpmax, &plotstart, &plotend);
159
160
161 cpgplot_identify(outf, obsexp, xypc, thresh, 0, len, shift,
162 ajStrGetPtr(bases), ajSeqGetNameC(seq),
163 minlen, minobsexp,
164 minpc, featout);
165
166 if(doplot)
167 {
168 ajGraphSetDatanameS(mult, ajSeqGetNameS(seq));
169 cpgplot_plotit(ajSeqGetNameC(seq), begin, len,
170 obsexp, xypc, thresh,
171 obsexpmax, doobsexp, docg, dopc, mult);
172 }
173
174 ajStrDel(&strand);
175 }
176
177 if(mult)
178 {
179 ajGraphicsClose();
180 ajGraphxyDel(&mult);
181 }
182
183 ajSeqDel(&seq);
184 ajStrDel(&substr);
185 ajStrDel(&bases);
186
187 ajFileClose(&outf);
188 ajSeqallDel(&seqall);
189 ajFeattabOutDel(&featout);
190
191 AJFREE(obsexp);
192 AJFREE(thresh);
193 AJFREE(xypc);
194
195 embExit();
196
197 return 0;
198 }
199
200
201
202
203 /* @funcstatic cpgplot_findbases **********************************************
204 **
205 ** Undocumented.
206 **
207 ** @param [r] substr [const AjPStr] Undocumented
208 ** @param [r] len [ajint] Undocumented
209 ** @param [r] window [ajint] Undocumented
210 ** @param [r] shift [ajint] Undocumented
211 ** @param [w] obsexp [float*] Undocumented
212 ** @param [w] xypc [float*] Undocumented
213 ** @param [r] bases [const AjPStr] Undocumented
214 ** @param [w] obsexpmax [float*] Undocumented
215 ** @param [w] plotstart [ajint*] Undocumented
216 ** @param [w] plotend [ajint*] Undocumented
217 ** @@
218 ******************************************************************************/
219
cpgplot_findbases(const AjPStr substr,ajint len,ajint window,ajint shift,float * obsexp,float * xypc,const AjPStr bases,float * obsexpmax,ajint * plotstart,ajint * plotend)220 static void cpgplot_findbases(const AjPStr substr, ajint len,
221 ajint window, ajint shift, float *obsexp,
222 float *xypc, const AjPStr bases,
223 float *obsexpmax,
224 ajint *plotstart, ajint *plotend)
225 {
226 float cxpy;
227 float cxf;
228 float cyf;
229 float windowf;
230
231
232 float obs;
233 float expect;
234 ajint i;
235 ajint j = 0;
236 ajint offset;
237
238 const char *p;
239 const char *q;
240
241 windowf = (float)window;
242 *obsexpmax = 0.0;
243 offset = window/2;
244 *plotstart = offset;
245 q = ajStrGetPtr(bases);
246
247 for(i=0; i<(len-window+1);i+=shift)
248 {
249 j = i+offset;
250 p = ajStrGetPtr(substr) + i;
251 cpgplot_countbases(p, q, window, &cxf, &cyf, &cxpy);
252
253 obs = cxpy;
254 expect = (cxf*cyf)/windowf;
255 if(!expect)
256 obsexp[j]=0.0;
257 else
258 {
259 obsexp[j] = obs/expect;
260 *obsexpmax = (*obsexpmax > obsexp[j]) ? *obsexpmax : obsexp[j];
261 }
262 xypc[j] = (cxf/windowf)*(float)100.0 + (cyf/windowf)*(float)100.0;
263 }
264
265 *plotend = j;
266
267 return;
268 }
269
270
271
272
273 /* @funcstatic cpgplot_countbases *********************************************
274 **
275 ** Undocumented.
276 **
277 ** @param [r] seq [const char*] Undocumented
278 ** @param [r] bases [const char*] Undocumented
279 ** @param [r] window [ajint] Undocumented
280 ** @param [w] cx [float*] Undocumented
281 ** @param [w] cy [float*] Undocumented
282 ** @param [w] cxpy [float*] Undocumented
283 ** @@
284 ******************************************************************************/
285
286
287
cpgplot_countbases(const char * seq,const char * bases,ajint window,float * cx,float * cy,float * cxpy)288 static void cpgplot_countbases(const char *seq, const char *bases,
289 ajint window,
290 float *cx, float *cy, float *cxpy)
291 {
292 ajint i;
293
294 ajint codex;
295 ajint codey;
296 ajint codea;
297 ajint codeb;
298
299 *cxpy = *cx = *cy = 0.;
300
301 codex = ajBaseAlphaToBin(bases[0]);
302 codey = ajBaseAlphaToBin(bases[1]);
303
304 codeb = ajBaseAlphaToBin(seq[0]);
305
306 for(i=0; i<window; ++i)
307 {
308 codea=codeb;
309 codeb=ajBaseAlphaToBin(seq[i+1]);
310
311 if(!(15-codea)) /* look for ambiguity code 'N' */
312 {
313 *cx = *cx + (float)0.25;
314 if(!(15-codeb))
315 *cxpy = *cxpy + (float)0.0625;
316 }
317 else
318 {
319 if(codea && !(codea & (15-codex)))
320 {
321 ++*cx;
322 if(codeb && !(codeb & (15-codey)))
323 ++*cxpy;
324 }
325 if(codea && !(codea & (15-codey)))
326 ++*cy;
327 }
328 }
329
330 return;
331 }
332
333
334
335
336 /* @funcstatic cpgplot_identify ***********************************************
337 **
338 ** This subroutine identifies the CpG line, identifying the possible
339 ** dinucleotide 'islands' in the sequence. These are defined as
340 ** base positions where, over an average of 10 windows, the calculated
341 ** % composition is over 50% and the calculated Obs/Exp ratio is over 0.6
342 ** and the conditions hold for a minimum of 200 bases.
343 **
344 ** @param [u] outf [AjPFile] Undocumented
345 ** @param [r] obsexp [const float*] Undocumented
346 ** @param [r] xypc [const float*] Undocumented
347 ** @param [w] thresh [AjBool*] Undocumented
348 ** @param [r] begin [ajint] Undocumented
349 ** @param [r] len [ajint] Undocumented
350 ** @param [r] shift [ajint] Undocumented
351 ** @param [r] bases [const char*] Undocumented
352 ** @param [r] name [const char*] Undocumented
353 ** @param [r] minlen [ajint] Undocumented
354 ** @param [r] minobsexp [float] Undocumented
355 ** @param [r] minpc [float] Undocumented
356 ** @param [u] featout [AjPFeattabOut] Undocumented
357 ** @@
358 ******************************************************************************/
359
cpgplot_identify(AjPFile outf,const float * obsexp,const float * xypc,AjBool * thresh,ajint begin,ajint len,ajint shift,const char * bases,const char * name,ajint minlen,float minobsexp,float minpc,AjPFeattabOut featout)360 static void cpgplot_identify(AjPFile outf,
361 const float *obsexp, const float *xypc,
362 AjBool *thresh, ajint begin, ajint len,
363 ajint shift, const char *bases, const char *name,
364 ajint minlen, float minobsexp, float minpc,
365 AjPFeattabOut featout)
366 {
367 static ajint avwindow = 10;
368 float avpc;
369 float avobsexp;
370 float sumpc;
371 float sumobsexp;
372
373 ajint i;
374 ajint pos;
375 ajint sumlen;
376 ajint first = 0;
377
378 for(i=0; i<len; ++i)
379 thresh[i] = ajFalse;
380
381 sumlen = 0;
382 for(pos=0,first=0;pos<(len-avwindow*shift);pos+=shift)
383 {
384 sumpc = sumobsexp = 0.0;
385 ajDebug("pos: %d max: %d\n", pos, pos+avwindow*shift);
386 for(i=pos;i<=(pos+avwindow*shift);++i)
387 {
388 ajDebug("obsexp[%d] %.2f xypc[%d] %.2f\n",
389 i, obsexp[i], i, xypc[i]);
390 sumpc += xypc[i];
391 sumobsexp += obsexp[i];
392 }
393
394 avpc = sumpc/(float)avwindow;
395 avobsexp = sumobsexp/(float)avwindow;
396 ajDebug("sumpc: %.2f sumobsexp: %.2f\n", sumpc, sumobsexp);
397 ajDebug(" avpc: %.2f avobsexp: %.2f\n", avpc, avobsexp);
398 if((avobsexp>minobsexp)&&(avpc>minpc))
399 {
400 if(!sumlen)
401 first = pos; /* start a new island */
402 sumlen += shift;
403 ajDebug(" ** hit first: %d sumlen: %d\n", first, sumlen);
404 }
405 else
406 {
407 if(sumlen >= minlen)
408 /* island long enough? */
409 for(i=first; i<=pos-shift;++i)
410 thresh[i] = ajTrue;
411
412 sumlen=0;
413 }
414 }
415
416 if(sumlen>=minlen)
417 for(i=first;i<len;++i)
418 thresh[i] = ajTrue;
419
420 cpgplot_reportislands(outf, thresh, bases, name, minobsexp, minpc,
421 minlen, begin, len);
422
423 cpgplot_dumpfeatout(featout,thresh, name, begin, len );
424
425 return;
426 }
427
428
429
430
431 /* @funcstatic cpgplot_reportislands ******************************************
432 **
433 ** Undocumented.
434 **
435 ** @param [u] outf [AjPFile] Undocumented
436 ** @param [r] thresh [const AjBool*] Undocumented
437 ** @param [r] bases [const char*] Undocumented
438 ** @param [r] name [const char*] Undocumented
439 ** @param [r] minobsexp [float] Undocumented
440 ** @param [r] minpc [float] Undocumented
441 ** @param [r] minlen [ajint] Undocumented
442 ** @param [r] begin [ajint] Undocumented
443 ** @param [r] len [ajint] Undocumented
444 ** @@
445 ******************************************************************************/
446
cpgplot_reportislands(AjPFile outf,const AjBool * thresh,const char * bases,const char * name,float minobsexp,float minpc,ajint minlen,ajint begin,ajint len)447 static void cpgplot_reportislands(AjPFile outf, const AjBool *thresh,
448 const char *bases, const char *name,
449 float minobsexp, float minpc,
450 ajint minlen, ajint begin, ajint len)
451 {
452 AjBool island;
453 ajint startpos = 0;
454 ajint endpos;
455 ajint slen;
456 ajint i;
457
458
459 ajFmtPrintF(outf,"\n\nCPGPLOT islands of unusual %s composition\n",
460 bases);
461 ajFmtPrintF(outf,"%s from %d to %d\n\n",name,begin+1,begin+len);
462 ajFmtPrintF(outf," Observed/Expected ratio > %.2f\n",minobsexp);
463 ajFmtPrintF(outf," Percent %c + Percent %c > %.2f\n",bases[0],
464 bases[1],minpc);
465 ajFmtPrintF(outf," Length > %d\n",minlen);
466
467 island = ajFalse;
468 for(i=0;i<len;++i)
469 {
470 if(island)
471 {
472 island = thresh[i];
473 if(!island)
474 {
475 slen = i - startpos;
476 endpos = i;
477 ajFmtPrintF(outf,"\n Length %d (%d..%d)\n",slen,
478 startpos+begin+1,endpos+begin);
479 }
480 }
481 else
482 {
483 island = thresh[i];
484 if(island)
485 startpos = i;
486 }
487 }
488
489 if(island)
490 {
491 slen = len-startpos+1;
492 endpos = len-1;
493 ajFmtPrintF(outf,"\n Length %d (%d..%d)\n",slen,
494 startpos+begin+1, endpos+begin);
495 }
496
497 return;
498 }
499
500
501
502
503 /* @funcstatic cpgplot_plotit *************************************************
504 **
505 ** Undocumented.
506 **
507 ** @param [r] seq [const char*] Undocumented
508 ** @param [r] begin [ajint] Undocumented
509 ** @param [r] len [ajint] Undocumented
510 ** @param [r] obsexp [const float*] Undocumented
511 ** @param [r] xypc [const float*] Undocumented
512 ** @param [r] thresh [const AjBool*] Undocumented
513 ** @param [r] obsexpmax [float] Undocumented
514 ** @param [r] doobsexp [AjBool] Undocumented
515 ** @param [r] docg [AjBool] Undocumented
516 ** @param [r] dopc [AjBool] Undocumented
517 ** @param [u] graphs [AjPGraph] Undocumented
518 ** @@
519 ******************************************************************************/
520
cpgplot_plotit(const char * seq,ajint begin,ajint len,const float * obsexp,const float * xypc,const AjBool * thresh,float obsexpmax,AjBool doobsexp,AjBool docg,AjBool dopc,AjPGraph graphs)521 static void cpgplot_plotit(const char *seq,
522 ajint begin, ajint len,
523 const float *obsexp, const float *xypc,
524 const AjBool *thresh,
525 float obsexpmax,
526 AjBool doobsexp, AjBool docg,
527 AjBool dopc, AjPGraph graphs)
528
529 {
530 AjPGraphdata tmGraph = NULL;
531 AjPGraphdata tmGraph2 = NULL;
532 AjPGraphdata tmGraph3 = NULL;
533 float *tmp = NULL;
534 ajint i;
535 float min = 0.;
536 float max = 0.;
537
538 ajint igraph=0;
539
540 if(doobsexp)
541 {
542 tmGraph2 = ajGraphdataNew();
543 ajGraphdataSetTitleC(tmGraph2,"Observed vs Expected");
544 ajGraphdataSetXlabelC(tmGraph2,"Base number");
545 ajGraphdataSetYlabelC(tmGraph2,"Obs/Exp");
546
547 min = 64000.;
548 max = -64000.;
549 for(i=0;i<len;++i)
550 {
551 min = (min<obsexp[i]) ? min : obsexp[i];
552 max = (max>obsexp[i]) ? max : obsexp[i];
553 }
554
555 ajGraphdataSetMinmax(tmGraph2,(float)begin,(float)begin+len-1,
556 min,max);
557 ajGraphdataSetTruescale(tmGraph2,(float)begin,(float)begin+len-1,
558 0.0,max);
559 ajGraphdataSetTypeC(tmGraph2,"Multi 2D Plot");
560
561 ajGraphxySetXstartF(graphs,(float)begin);
562 ajGraphxySetXendF(graphs,(float)(begin+len-1));
563
564 ajGraphxySetYstartF(graphs,0.0);
565 ajGraphxySetYendF(graphs,obsexpmax);
566 ajGraphxySetXrangeII(graphs,begin,begin+len-1);
567 ajGraphxySetYrangeII(graphs,0,(ajint)(obsexpmax+1.0));
568
569 ajGraphdataCalcXY(tmGraph2,len,(float)begin,1.0,obsexp);
570 ajGraphDataReplaceI(graphs,tmGraph2,igraph++);
571 tmGraph2 = NULL;
572 }
573
574 if(dopc)
575 {
576 tmGraph3 = ajGraphdataNew();
577 ajGraphdataSetTitleC(tmGraph3,"Percentage");
578 ajGraphdataSetXlabelC(tmGraph3,"Base number");
579 ajGraphdataSetYlabelC(tmGraph3,"Percentage");
580
581 min = 64000.;
582 max = -64000.;
583 for(i=0;i<len;++i)
584 {
585 min = (min<xypc[i]) ? min : xypc[i];
586 max = (max>xypc[i]) ? max : xypc[i];
587 }
588
589 ajGraphdataSetMinmax(tmGraph3,(float)begin,(float)begin+len-1,
590 min,max);
591 ajGraphdataSetTruescale(tmGraph3,(float)begin,(float)begin+len-1,
592 0.0,max);
593 ajGraphdataSetTypeC(tmGraph3,"Multi 2D Plot");
594
595 ajGraphxySetXstartF(graphs,(float)begin);
596 ajGraphxySetXendF(graphs,(float)(begin+len-1));
597 ajGraphxySetYstartF(graphs,0);
598 ajGraphxySetYendF(graphs,100);
599 ajGraphxySetXrangeII(graphs,begin,begin+len-1);
600 ajGraphxySetYrangeII(graphs,0,100);
601
602 ajGraphdataCalcXY(tmGraph3,len,(float)begin,1.0,xypc);
603 ajGraphDataReplaceI(graphs,tmGraph3,igraph++);
604 tmGraph3 = NULL;
605 }
606
607 if(docg)
608 {
609 AJCNEW (tmp, len);
610 for(i=0;i<len;++i)
611 {
612 if(thresh[i])
613 tmp[i] = 1.0;
614 else
615 tmp[i] = 0.0;
616 }
617
618 tmGraph = ajGraphdataNew();
619 ajGraphdataSetTitleC(tmGraph,"Putative Islands");
620 ajGraphdataSetXlabelC(tmGraph,"Base Number");
621 ajGraphdataSetYlabelC(tmGraph,"Threshold");
622
623 ajGraphdataSetMinmax(tmGraph,(float)begin,(float)begin+len-1,
624 0.,1.);
625 ajGraphdataSetTypeC(tmGraph,"Multi 2D Plot");
626
627
628 ajGraphxySetXstartF(graphs,(float)begin);
629 ajGraphxySetXendF(graphs,(float)(begin+len-1));
630 ajGraphxySetYstartF(graphs,0);
631 ajGraphxySetYendF(graphs,2);
632 ajGraphxySetXrangeII(graphs,begin,begin+len-1);
633 ajGraphxySetYrangeII(graphs,0,2);
634 ajGraphdataSetMinmax(tmGraph,(float)begin,(float)(begin+len-1),
635 (float) 0.0, (float) 1.2);
636 ajGraphdataSetTruescale(tmGraph,(float)begin,(float)(begin+len-1),
637 (float) 0.0, (float) 1.0);
638
639 ajGraphdataCalcXY(tmGraph,len,(float)begin,1.0,tmp);
640 ajGraphDataReplaceI(graphs,tmGraph,igraph++);
641 tmGraph = NULL;
642 }
643
644
645 if(docg || dopc || doobsexp)
646 {
647 ajGraphSetTitleC(graphs,seq);
648 ajGraphxySetflagOverlay(graphs,ajFalse);
649 ajGraphxyDisplay(graphs, AJFALSE);
650 }
651
652
653 if(docg)
654 AJFREE (tmp);
655
656 return;
657 }
658
659
660
661
662 /* @funcstatic cpgplot_dumpfeatout ********************************************
663 **
664 ** Undocumented.
665 **
666 ** @param [u] featout [AjPFeattabOut] Undocumented
667 ** @param [r] thresh [const AjBool*] Undocumented
668 ** @param [r] seqname [const char*] Undocumented
669 ** @param [r] begin [ajint] Undocumented
670 ** @param [r] len [ajint] Undocumented
671 ** @@
672 ******************************************************************************/
673
cpgplot_dumpfeatout(AjPFeattabOut featout,const AjBool * thresh,const char * seqname,ajint begin,ajint len)674 static void cpgplot_dumpfeatout(AjPFeattabOut featout, const AjBool *thresh,
675 const char *seqname, ajint begin, ajint len)
676 {
677 AjBool island;
678 ajint startpos = 0;
679 ajint endpos;
680 ajint i;
681 AjPFeattable feattable;
682 AjPStr name = NULL;
683 AjPStr source = NULL;
684 AjPStr type = NULL;
685 char strand = '+';
686 ajint frame = 0;
687 AjPFeature feature;
688 float score = 0.0;
689
690 ajStrAssignC(&name,seqname);
691
692 feattable = ajFeattableNewDna(name);
693
694 ajStrAssignC(&source,"cpgplot");
695 ajStrAssignC(&type,"misc_feature");
696
697
698 island = ajFalse;
699 for(i=0;i<len;++i)
700 {
701 if(island)
702 {
703 island = thresh[i];
704 if(!island)
705 {
706 endpos = i;
707 feature = ajFeatNew(feattable, source, type,
708 startpos+begin+1,endpos+begin,
709 score, strand, frame) ;
710 if(!feature)
711 ajDebug("Error feature not added to feature table");
712 }
713 }
714 else
715 {
716 island = thresh[i];
717 if(island)
718 startpos = i;
719 }
720 }
721
722 if(island)
723 {
724 endpos = len-1;
725 feature = ajFeatNew(feattable, source, type,
726 startpos+begin+1,endpos+begin,
727 score, strand, frame);
728 }
729 ajFeatSortByStart(feattable);
730 ajFeattableWrite (featout, feattable);
731 ajFeattableDel(&feattable);
732
733 ajStrDel(&source);
734 ajStrDel(&type);
735 ajStrDel(&name);
736
737 return;
738 }
739