1 /* @source recoder
2 **
3 ** Find restriction sites in a nucleotide sequence and remove
4 ** them whilst maintaining the same translation.
5 **
6 **
7 ** This program uses a similar output style to that of
8 ** "silent". Also, the enzyme reading function was taken
9 ** straight out of "silent".
10 **
11 ** @author Copyright (C) Tim Carver (tcarver@hgmp.mrc.ac.uk)
12 ** @@
13 **
14 ** This program is free software; you can redistribute it and/or
15 ** modify it under the terms of the GNU General Public License
16 ** as published by the Free Software Foundation; either version 2
17 ** of the License, or (at your option) any later version.
18 **
19 ** This program is distributed in the hope that it will be useful,
20 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
21 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 ** GNU General Public License for more details.
23 **
24 ** You should have received a copy of the GNU General Public License
25 ** along with this program; if not, write to the Free Software
26 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
27 ******************************************************************************/
28
29 #include "emboss.h"
30 #include <limits.h>
31
32 #define IUBFILE "Ebases.iub"
33
34
35
36
37 /* @datastatic PRinfo *******************************************************
38 **
39 ** recoder internals for RE information
40 **
41 ** @alias SRinfo
42 ** @alias ORinfo
43 **
44 ** @attr code [AjPStr] Undocumented
45 ** @attr site [AjPStr] Undocumented
46 ** @attr revsite [AjPStr] Undocumented
47 ** @attr ncuts [ajint] Undocumented
48 ** @attr cut1 [ajint] Undocumented
49 ** @attr cut2 [ajint] Undocumented
50 ** @attr cut3 [ajint] Undocumented
51 ** @attr cut4 [ajint] Undocumented
52 ** @attr Padding [char[4]] Padding to alignment boundary
53 ******************************************************************************/
54
55 typedef struct SRinfo
56 {
57 AjPStr code;
58 AjPStr site;
59 AjPStr revsite;
60 ajint ncuts;
61 ajint cut1;
62 ajint cut2;
63 ajint cut3;
64 ajint cut4;
65 char Padding[4];
66 } ORinfo;
67 #define PRinfo ORinfo*
68
69
70
71
72 /* @datastatic PMutant *********************************************************
73 **
74 ** recoder internals for mutation sites
75 **
76 ** @alias SMutant
77 ** @alias OMutant
78 **
79 ** @attr code [AjPStr] Undocumented
80 ** @attr site [AjPStr] Undocumented
81 ** @attr match [ajint] Undocumented
82 ** @attr base [ajint] Undocumented
83 ** @attr seqaa [AjPStr] Undocumented
84 ** @attr reaa [AjPStr] Undocumented
85 ** @attr obase [char] Undocumented
86 ** @attr nbase [char] Undocumented
87 ** @attr Padding [char[6]] Padding to alignment boudnary
88 ******************************************************************************/
89
90 typedef struct SMutant
91 {
92 AjPStr code;
93 AjPStr site;
94 ajint match;
95 ajint base;
96 AjPStr seqaa;
97 AjPStr reaa;
98 char obase;
99 char nbase;
100 char Padding[6];
101 } OMutant;
102 #define PMutant OMutant*
103
104
105 static AjPTrn recoderTable = NULL; /* translation table object */
106
107
108 static ajint recoder_readRE(AjPList *relist, const AjPStr enzymes);
109 static AjPList recoder_rematch(const AjPStr sstr, AjPList ressite,
110 AjPStr* tailstr,
111 const AjPStr sname, ajint RStotal,
112 ajint begin, ajint end,
113 AjBool tshow);
114 static AjPList recoder_checkTrans(const AjPStr seq,const EmbPMatMatch match,
115 const PRinfo rlp, ajint begin,
116 ajint pos, AjBool rev,
117 AjBool* empty);
118 static AjBool recoder_checkPat(const EmbPMatMatch match,
119 const PRinfo rlp,
120 ajint end);
121 static ajint recoder_changebase(char pbase, char* tbase);
122 static void recoder_mutFree(PMutant* mut);
123 static ajint recoder_basecompare(const void *a, const void *b);
124
125 static void recoder_fmt_seq(const char* title, const AjPStr seq,
126 AjPStr* tailstr,
127 ajint start, AjBool num);
128 static void recoder_fmt_muts(AjPList muts, AjPFeattable feat);
129
130
131
132
133 /* @prog recoder **************************************************************
134 **
135 ** Remove restriction sites but maintain the same translation
136 **
137 ******************************************************************************/
138
main(int argc,char ** argv)139 int main(int argc, char **argv)
140 {
141 AjPSeq seq = NULL;
142 AjPReport report = NULL;
143 AjPFeattable feat=NULL;
144
145 AjPStr sstr = NULL;
146 const AjPStr sname = NULL; /* seq name */
147 AjPStr enzymes = NULL; /* string for RE selection */
148
149 ajint RStotal;
150 ajint begin; /* specified start by user */
151 ajint end; /* specified end by user */
152 ajint start;
153 AjBool sshow;
154 AjBool tshow;
155
156 AjPList relist = NULL;
157 AjPList muts;
158 PRinfo re;
159 AjPStr tailstr = NULL;
160
161 embInit("recoder", argc, argv);
162
163 seq = ajAcdGetSeq("sequence"); /* sequence to investigate */
164 enzymes = ajAcdGetString("enzymes"); /* enzyme list */
165 sshow = ajAcdGetBoolean("sshow"); /* display seq */
166 tshow = ajAcdGetBoolean("tshow"); /* display translated seq */
167 report = ajAcdGetReport ("outfile"); /* report filename */
168
169
170 RStotal=recoder_readRE(&relist,enzymes); /* read in RE info */
171
172 ajDebug("readRE read %d listlen:%Lu\n", RStotal, ajListGetLength(relist));
173
174 begin = ajSeqGetBegin(seq); /* seq start posn, or 1 */
175 end = ajSeqGetEnd(seq); /* seq end posn, or seq len */
176
177 /* --begin and --end to convert counting from 0-N, not 1-N */
178 ajStrAssignSubC(&sstr,ajSeqGetSeqC(seq),--begin,--end);
179 ajStrFmtUpper(&sstr);
180
181 sname = ajSeqGetNameS(seq);
182
183 start = begin+1;
184
185 feat = ajFeattableNewDna(ajSeqGetNameS(seq));
186
187 if(sshow)
188 {
189 recoder_fmt_seq("SEQUENCE",
190 sstr,&tailstr,start,ajTrue);
191 }
192
193 /******* get de-restriction site *******/
194 /******* forward strand *******/
195 muts = recoder_rematch(sstr,relist,&tailstr,sname,RStotal,
196 begin,end,tshow);
197
198
199 ajReportSetHeaderC(report,
200 "KEY:\n"
201 "EnzymeName: Enzyme name\n"
202 "RS-Pattern: Restriction enzyme recognition site "
203 "pattern\n"
204 "Base-Posn: Position of base to be mutated\n"
205 "AAs: Amino acid. Original sequence(.)After mutation\n"
206 "Mutation: The base mutation to perform\n\n"
207 "Creating silent mutations");
208 recoder_fmt_muts(muts,feat);
209
210 ajReportSetTailS(report, tailstr);
211 ajReportSetStatistics(report, 1, ajSeqGetLenTrimmed(seq));
212 (void) ajReportWrite (report,feat,seq);
213 ajFeattableDel(&feat);
214
215 while(ajListPop(relist,(void **)&re))
216 {
217 ajStrDel(&re->code);
218 ajStrDel(&re->site);
219 ajStrDel(&re->revsite);
220 AJFREE(re);
221 }
222 ajListFree(&relist);
223
224 ajSeqDel(&seq);
225 ajStrDel(&enzymes);
226 ajStrDel(&sstr);
227 ajStrDel(&tailstr);
228
229 ajListFree(&muts);
230
231 ajReportClose(report);
232 ajReportDel(&report);
233 ajTrnDel(&recoderTable);
234
235 embExit();
236
237 return 0;
238 }
239
240
241
242
243 /* @funcstatic recoder_rematch ************************************************
244 **
245 ** Looks for RE matches and test new bases to find those that
246 ** give the same translation.
247 **
248 ** @param [r] sstr [const AjPStr] Search sequence as a string
249 ** @param [u] relist [AjPList] Regular expression list
250 ** @param [w] tailstr [AjPStr*] Report tail as a string
251 ** @param [r] sname [const AjPStr] Sequence name
252 ** @param [r] RStotal [ajint] Restriction sites
253 ** @param [r] begin [ajint] Start position
254 ** @param [r] end [ajint] End position
255 ** @param [r] tshow [AjBool] Show translation
256 ** @return [AjPList] Mutant list of those RS sites removed but
257 ** maintaining same translation.
258 ******************************************************************************/
259
recoder_rematch(const AjPStr sstr,AjPList relist,AjPStr * tailstr,const AjPStr sname,ajint RStotal,ajint begin,ajint end,AjBool tshow)260 static AjPList recoder_rematch(const AjPStr sstr, AjPList relist,
261 AjPStr* tailstr,
262 const AjPStr sname, ajint RStotal,
263 ajint begin, ajint end,
264 AjBool tshow)
265 {
266 AjPList res;
267 AjPList results;
268 AjPStr str; /* holds RS patterns */
269 AjPStr tstr;
270 AjPStr pep = NULL; /* string to hold protein */
271
272 AjBool dummy=ajFalse; /* need a bool for ajPatClassify */
273 AjBool empty;
274
275 ajint mm; /* no. of mismatches */
276 ajint pats; /* no. of pattern matches */
277 ajint patlen;
278 ajint pos;
279 ajint aw;
280 ajint start;
281
282 AjPList patlist = NULL; /* list for pattern matches of.. */
283 EmbPMatMatch match; /* ..AjMatMatch structures*/
284 PRinfo rlp = NULL;
285
286 AjBool rev = ajFalse;
287
288 str = ajStrNew();
289 tstr = ajStrNew();
290 pep = ajStrNew();
291 if(!recoderTable)
292 recoderTable = ajTrnNewI(0); /* 0 for std DNA (see fuzztran) */
293
294 results = ajListNew();
295 start = 1;
296
297 ajTrnSeqFrameS(recoderTable,sstr,1,&pep); /* frame 1 */
298 if(tshow)
299 {
300 recoder_fmt_seq("TRANSLATED SEQUENCE",
301 pep,tailstr,start,ajFalse);
302 }
303
304 for(aw=0;aw<RStotal;aw++) /* read in RE patterns */
305 {
306 /* pop off first RE */
307 (void) ajListPop(relist,(void **)&rlp);
308 /* ignore unknown cut sites & zero cutters */
309 if(*ajStrGetPtr(rlp->site)=='?'||!rlp->ncuts)
310 {
311 ajListPushAppend(relist,(void*)rlp);
312 continue;
313 }
314
315 ajStrFmtUpper(&rlp->site); /* RS to upper case */
316 ajStrFmtUpper(&rlp->revsite);
317
318 ajStrAssignS(&str,rlp->site); /* str holds RS pat */
319
320 rev = ajFalse;
321
322 while(str)
323 {
324 patlist = ajListNew(); /* list for matches */
325
326 /* convert our RS pattern to a reg expression */
327 /* "0" means DNA & so does any ambig conversions */
328
329 if(!embPatClassify(str,&str,&dummy,&dummy,&dummy,
330 &dummy,&dummy,&dummy,0)) continue;
331
332 /* find pattern matches in seq with NO mismatches */
333 mm = 0;
334 pats = embPatBruteForce(sstr,str,ajFalse,ajFalse,
335 patlist,begin+1,mm,sname);
336
337 if(pats)
338 {
339 while(ajListPop(patlist,(void**)&match))
340 {
341 patlen = match->len;
342
343 if(recoder_checkPat(match,rlp,
344 end))
345 {
346 for(pos=0;pos<patlen;pos++)
347 {
348 res = recoder_checkTrans(sstr,match,rlp,begin,
349 pos,rev, &empty);
350 if(empty)
351 ajListFree(&res);
352 else
353 ajListPushlist(results,&res);
354 }
355 }
356 embMatMatchDel(&match);
357 }
358 }
359
360 ajListFree(&patlist);
361
362 if(!rev && !ajStrMatchS(rlp->site, rlp->revsite))
363 {
364 ajStrAssignS(&str,rlp->revsite); /* str holds RS pat */
365 rev = ajTrue;
366 }
367 else
368 ajStrDel(&str);
369 }
370 /* push RE info back to top of the list */
371 ajListPushAppend(relist,(void *)rlp);
372 }
373
374 ajStrDel(&str);
375 ajStrDel(&tstr);
376 ajStrDel(&pep);
377
378 return results;
379 }
380
381
382
383
384 /* @funcstatic recoder_readRE *************************************************
385 **
386 ** Read in RE information from REBASE file.
387 **
388 ** @param [w] relist [AjPList*] returns restriction site information as a list
389 ** @param [r] enzymes [const AjPStr] Selected enzymes to read
390 ** @return [ajint] Number of restriction sites in list
391 **
392 ******************************************************************************/
recoder_readRE(AjPList * relist,const AjPStr enzymes)393 static ajint recoder_readRE(AjPList *relist, const AjPStr enzymes)
394 {
395 EmbPPatRestrict rptr = NULL; /* store RE info */
396 AjPFile fin = NULL; /* file pointer to RE file data */
397 AjPStr refilename = NULL; /* .. & string for the filename */
398 register ajint RStotal = 0; /* counts no of RE */
399 PRinfo rinfo = NULL;
400 AjBool isall = ajFalse;
401 ajint ne = 0;
402 ajint i;
403 AjPStr *ea = NULL;
404
405 refilename = ajStrNewC("REBASE/embossre.enz");
406
407 rptr = embPatRestrictNew(); /* allocate a restrict struc */
408
409 *relist = ajListNew(); /* list the RS code and info */
410
411 fin = ajDatafileNewInNameS(refilename);
412 if(!fin)
413 ajFatal("Aborting...restriction file not found");
414
415 if(!enzymes) /* parse user-selected enzyme list */
416 isall = ajTrue;
417 else
418 {
419 ne=ajArrCommaList(enzymes,&ea); /* no. of RE's specified */
420 /* ea points to enzyme list */
421 for(i=0;i<ne;++i)
422 ajStrRemoveWhite(&ea[i]); /* remove all whitespace */
423
424 if(ajStrMatchCaseC(ea[0],"all"))
425 isall = ajTrue;
426 else
427 isall = ajFalse;
428 }
429 ajDebug("recoder_readRE enzymes: '%S' isall:%B\n", enzymes, isall);
430
431
432 /* read RE data into AjPRestrict obj */
433
434 while(!ajFileIsEof(fin))
435 {
436 if(!embPatRestrictReadEntry(rptr,fin))
437 continue;
438
439 if(!isall) /* only select enzymes on command line */
440 {
441 for(i=0;i<ne;++i)
442 if(ajStrMatchCaseS(ea[i],rptr->cod))
443 break;
444
445 if(i==ne)
446 continue;
447 }
448
449
450 AJNEW(rinfo);
451 rinfo->code = ajStrNewS(rptr->cod);
452 rinfo->site = ajStrNewS(rptr->pat);
453 rinfo->revsite = ajStrNewS(rptr->pat);
454 ajSeqstrReverse(&rinfo->revsite);
455 rinfo->ncuts = rptr->ncuts;
456 rinfo->cut1 = rptr->cut1;
457 rinfo->cut2 = rptr->cut2;
458 rinfo->cut3 = rptr->cut3;
459 rinfo->cut4 = rptr->cut4;
460 ajListPush(*relist,(void *)rinfo);
461 RStotal++;
462 ajDebug("Creating RStotal:%d listsize:%Lu '%S' '%S'\n",
463 RStotal, ajListGetLength(*relist), rptr->cod, rptr->pat);
464 }
465
466 for(i=0;i<ne;++i)
467 ajStrDel(&ea[i]);
468 AJFREE(ea);
469 embPatRestrictDel(&rptr);
470 ajFileClose(&fin);
471 ajStrDel(&refilename);
472
473 return RStotal;
474 }
475
476
477
478
479 /* @funcstatic recoder_checkPat ***********************************************
480 **
481 ** Checks whether the RS pattern falls within the sequence string
482 **
483 ** @param [r] match [const EmbPMatMatch] Match data
484 ** @param [r] rlp [const PRinfo] Restriction site info
485 ** @param [r] end [ajint] End position
486 ** @return [AjBool] ajTrue if the pattern is found
487 **
488 ******************************************************************************/
recoder_checkPat(const EmbPMatMatch match,const PRinfo rlp,ajint end)489 static AjBool recoder_checkPat(const EmbPMatMatch match,
490 const PRinfo rlp, ajint end)
491 {
492 ajint mpos;
493
494 ajint min = INT_MAX; /* reverse sense intentional! */
495 ajint max = -INT_MAX;
496
497
498 mpos = match->start; /* start posn of match in seq */
499
500 if(rlp->ncuts==4) /* test if cut site is within seq */
501 {
502 min = AJMIN(rlp->cut1,rlp->cut2);
503 max = AJMAX(rlp->cut3,rlp->cut4);
504 }
505 else if(rlp->ncuts==2)
506 {
507 min = AJMIN(rlp->cut1,rlp->cut2);
508 max = AJMAX(rlp->cut1,rlp->cut2);
509 }
510 else
511 {
512 ajWarn("Possibly corrupt RE file");
513 return ajFalse;
514 }
515
516 if(mpos+min<0||mpos+max>end+1)
517 return ajFalse; /* cut not in seq */
518
519 return ajTrue;
520 }
521
522
523
524
525 /* @funcstatic recoder_checkTrans *********************************************
526 **
527 ** Identify mutations at a site in the RS pattern that result in the
528 ** same translation.
529 **
530 ** @param [r] dna [const AjPStr] Sequence as a string
531 ** @param [r] match [const EmbPMatMatch] Match data
532 ** @param [r] rlp [const PRinfo] Restriction site info
533 ** @param [r] begin [ajint] Start position
534 ** @param [r] pos [ajint] Base position in site
535 ** @param [r] rev [AjBool] Use reverse site
536 ** @param [w] empty [AjBool*] ajTrue if the list is empty
537 ** @return [AjPList] List of de-restricted sites which maintain same
538 ** translation.
539 **
540 ******************************************************************************/
541
recoder_checkTrans(const AjPStr dna,const EmbPMatMatch match,const PRinfo rlp,ajint begin,ajint pos,AjBool rev,AjBool * empty)542 static AjPList recoder_checkTrans(const AjPStr dna, const EmbPMatMatch match,
543 const PRinfo rlp, ajint begin,
544 ajint pos, AjBool rev, AjBool* empty)
545 {
546 char *pseq;
547 char *pseqm;
548 const char *prs;
549 char *s;
550
551 PMutant tresult;
552 AjPList res;
553
554 ajint mpos;
555 ajint i;
556 AjPStr s1 = NULL;
557 AjPStr s2 = NULL;
558 char base;
559 char pbase;
560 char tbase[4];
561
562 ajint x;
563 ajint nb;
564
565 AjPStr seq = NULL;
566
567 seq = ajStrNewS(dna);
568
569 *empty = ajTrue;
570 mpos = match->start; /* start posn of match in seq */
571
572 pseq = ajStrGetuniquePtr(&seq); /* pointer to start of seq */
573 pseqm = pseq+mpos-(begin+1); /* pointer to start of match in seq */
574
575 if(!rev)
576 prs = ajStrGetPtr(rlp->site); /* pointer to start of RS pattern */
577 else
578 prs = ajStrGetPtr(rlp->revsite);
579
580 base = pseqm[pos]; /* store orig seq base */
581 pbase = prs[pos];
582
583 /* use IUB codes to get other bases */
584 nb = recoder_changebase(pbase,&tbase[0]);
585
586 x = mpos+pos-(begin+1);
587
588 s = pseq+x-x%3;
589
590 if(!recoderTable)
591 recoderTable = ajTrnNewI(0);
592 s1 = ajStrNewK(ajTrnCodonC(recoderTable,s));
593
594 res=ajListNew();
595
596 for(i=0;i<nb;i++) /* try out other bases */
597 {
598 pseq[x] = tbase[i];
599 s2 = ajStrNewK(ajTrnCodonC(recoderTable,s));
600
601 if(ajStrMatchS(s1,s2))
602 {
603 /* if same translation */
604 AJNEW(tresult);
605 tresult->obase = base;
606 tresult->nbase = tbase[i];
607 tresult->code = ajStrNewS(rlp->code);
608 if(rev)
609 tresult->site = ajStrNewS(rlp->revsite);
610 else
611 tresult->site = ajStrNewS(rlp->site);
612 tresult->seqaa = ajStrNewS(s1);
613 tresult->reaa = ajStrNewS(s2);
614
615 tresult->match = mpos;
616 tresult->base = mpos+pos;
617
618 ajListPushAppend(res,(void *)tresult);
619 *empty = ajFalse;
620 }
621 ajStrDel(&s2);
622 }
623
624 pseq[x] = base; /* resubstitute orig base */
625
626 ajStrDel(&s1);
627 ajStrDel(&seq);
628
629 return res;
630 }
631
632
633
634
635 /* @funcstatic recoder_changebase *********************************************
636 **
637 ** Use IUB code to return alternative nucleotides to that provided
638 ** which result in the same translation.
639 **
640 ** @param [r] pbase [char] Base
641 ** @param [w] tbase [char*] C string with alternative bases
642 ** @return [ajint] number of bases stored in tbase
643 **
644 ******************************************************************************/
645
recoder_changebase(char pbase,char * tbase)646 static ajint recoder_changebase(char pbase, char* tbase)
647 {
648 ajint setBase[] =
649 {
650 1,1,1,1
651 };
652 AjIStr splits = NULL;
653 const AjPStr bt = NULL;
654 char bs;
655 ajint i;
656 ajint nb;
657
658 bt = ajBaseGetCodes((ajint)pbase);
659 splits = ajStrIterNew(bt);
660
661 while(!ajStrIterDone(splits))
662 {
663 bs = ajStrIterGetK(splits);
664
665 if( ajBaseAlphaToBin(bs) & ajBaseAlphaToBin('G') )
666 setBase[0] = 0;
667 if( ajBaseAlphaToBin(bs) & ajBaseAlphaToBin('A') )
668 setBase[1] = 0;
669 if( ajBaseAlphaToBin(bs) & ajBaseAlphaToBin('T') )
670 setBase[2] = 0;
671 if( ajBaseAlphaToBin(bs) & ajBaseAlphaToBin('C') )
672 setBase[3] = 0;
673
674 ajStrIterNext(splits);
675 }
676
677 ajStrIterDel(&splits);
678
679 nb = 0;
680 for(i=0;i<4;i++)
681 {
682 if( setBase[i] == 1 )
683 {
684 if(i==0)
685 tbase[nb] = 'G';
686 else if(i==1)
687 tbase[nb] = 'A';
688 else if(i==2)
689 tbase[nb] = 'T';
690 else if(i==3)
691 tbase[nb] = 'C';
692 nb++;
693 }
694 }
695
696 return nb;
697 }
698
699
700
701
702 /* @funcstatic recoder_fmt_seq ************************************************
703 **
704 ** Write sequence to the output file.
705 **
706 ** @param [r] title [const char*] Title for sequence report
707 ** @param [r] seq [const AjPStr] Sequence as a string
708 ** @param [w] tailstr [AjPStr*] Report tail as a string
709 ** @param [r] start [ajint] Start position
710 ** @param [r] num [AjBool] Numbered sequence
711 ** @return [void]
712 **
713 ******************************************************************************/
714
recoder_fmt_seq(const char * title,const AjPStr seq,AjPStr * tailstr,ajint start,AjBool num)715 static void recoder_fmt_seq(const char* title, const AjPStr seq,
716 AjPStr* tailstr,
717 ajint start, AjBool num)
718 {
719 const char *p;
720 ajint m;
721 ajint i;
722 ajint tlen;
723
724 ajFmtPrintAppS(tailstr,"%s:\n",title);
725 if(num)
726 {
727 p=ajStrGetPtr(seq);
728 ajFmtPrintAppS(tailstr,"%-7d",start);
729 tlen=ajStrGetLen(seq);
730 for(i=0; i<tlen ; i++)
731 {
732 ajFmtPrintAppS(tailstr,"%c",p[i]);
733 m=i+1;
734 if(m%10==0)
735 ajFmtPrintAppS(tailstr," ");
736 if(m%60==0 && m<tlen)
737 ajFmtPrintAppS(tailstr,"\n%-7d",(start+m));
738 }
739 }
740 else
741 {
742 p=ajStrGetPtr(seq);
743 tlen=ajStrGetLen(seq);
744 for(i=0; i<tlen ; i++)
745 {
746 ajFmtPrintAppS(tailstr,"%c",p[i]);
747 m=i+1;
748 if(m%10==0)
749 ajFmtPrintAppS(tailstr," ");
750 if(m%60==0)
751 ajFmtPrintAppS(tailstr,"\n");
752 }
753 }
754
755 ajFmtPrintAppS(tailstr,"\n\n");
756
757 return;
758 }
759
760
761
762
763 /* @funcstatic recoder_fmt_muts ***********************************************
764 **
765 ** Write de-restricted sites to feature table
766 **
767 ** @param [u] muts [AjPList] List of derestricted sites
768 ** @param [u] feat [AjPFeattable] Feature table object
769 ** @return [void]
770 ******************************************************************************/
771
recoder_fmt_muts(AjPList muts,AjPFeattable feat)772 static void recoder_fmt_muts(AjPList muts, AjPFeattable feat)
773 {
774 PMutant res;
775 AjPFeature sf = NULL;
776 AjPStr tmpFeatStr = NULL;
777
778 ajListSort(muts, &recoder_basecompare);
779
780 while(ajListPop(muts,(void **)&res))
781 {
782 sf = ajFeatNewII(feat,
783 res->match, res->match+ajStrGetLen(res->site)-1);
784
785 ajFmtPrintS(&tmpFeatStr, "*enzyme %S", res->code);
786 ajFeatTagAddSS (sf, NULL, tmpFeatStr);
787 ajFmtPrintS(&tmpFeatStr, "*rspattern %S", res->site);
788 ajFeatTagAddSS(sf, NULL, tmpFeatStr);
789 ajFmtPrintS(&tmpFeatStr, "*baseposn %d", res->base);
790 ajFeatTagAddSS(sf, NULL, tmpFeatStr);
791 ajFmtPrintS(&tmpFeatStr, "*aa %S.%S", res->seqaa, res->reaa);
792 ajFeatTagAddSS(sf, NULL, tmpFeatStr);
793 ajFmtPrintS(&tmpFeatStr, "*mutation %c->%c", res->obase,res->nbase);
794 ajFeatTagAddSS(sf, NULL, tmpFeatStr);
795
796 recoder_mutFree(&res);
797 }
798
799 ajStrDel(&tmpFeatStr);
800 return;
801 }
802
803
804
805
806 /* @funcstatic recoder_basecompare ********************************************
807 **
808 ** Compare 2 base positions in the nucleotide sequence
809 **
810 ** @param [r] a [const void *] First base
811 ** @param [r] b [const void *] Second base
812 ** @return [ajint] Comparison result
813 ******************************************************************************/
814
recoder_basecompare(const void * a,const void * b)815 static ajint recoder_basecompare(const void *a, const void *b)
816 {
817 return((*(PMutant const *)a)->base)-((*(PMutant const *)b)->base);
818 }
819
820
821
822
823 /* @funcstatic recoder_mutFree ************************************************
824 **
825 ** Free allocated memory for mutant structure
826 **
827 ** @param [d] mut [PMutant*] Mutant structure to be deleted
828 ** @return [void]
829 ******************************************************************************/
830
recoder_mutFree(PMutant * mut)831 static void recoder_mutFree(PMutant* mut)
832 {
833 ajStrDel(&(*mut)->code);
834 ajStrDel(&(*mut)->site);
835 ajStrDel(&(*mut)->seqaa);
836 ajStrDel(&(*mut)->reaa);
837 AJFREE(*mut);
838
839 return;
840 }
841