diff -Naur emboss/acd/silent.acd emboss_silent/acd/silent.acd --- emboss/acd/silent.acd 2003-07-30 09:25:56.000000000 +0000 +++ emboss_silent/acd/silent.acd 2003-08-05 07:51:18.000000000 +0000 @@ -32,23 +32,40 @@ type: "page" ] - boolean: sshow [ + boolean: complement [ default: "N" - information: "Display untranslated sequence" + information: "Display complementary strand" ] - boolean: tshow [ + boolean: sshow [ default: "N" - information: "Display translated sequence" + information: "Display untranslated sequence" ] + # boolean: tshow [ + # default: "N" + # information: "Display translated sequence" + # ] + boolean: allmut [ default: "N" information: "Display all mutations" ] - outfile: outf [ + report: outrep [ parameter: "Y" + information: "Reports file name" + rformat: "table" + rscoreshow: "N" + multiple: "N" + taglist: "str:enzyme=Enzyme str:site=Site str:origaa=OrigAA + str:mutaa=MutAA int:base=Base str:orign=OrigN + str:mutn=MutN" + ] + + seqout: outseq [ + extension: "seq" + information: "Sequences file name" ] endsection: output diff -Naur emboss/silent.c emboss_silent/silent.c --- emboss/silent.c 2003-07-30 09:25:56.000000000 +0000 +++ emboss_silent/silent.c 2003-08-05 07:51:52.000000000 +0000 @@ -52,24 +52,28 @@ /* Prototypes */ -static AjPList silent_mismatch(AjPStr sstr, AjPList ressite, AjPFile outf, +// static AjPList silent_mismatch(AjPStr sstr, AjPList ressite, AjPFile outf, +static AjPList silent_mismatch(AjPStr sstr, AjPList ressite, AjPSeqout outseq, AjPStr sname, ajint RStotal, ajint begin, ajint radj, AjBool rev, ajint end, AjBool tshow); static ajint silent_restr_read(AjPList *relist, AjPStr enzymes); -static AjBool silent_checktrans(AjPStr seq,AjPFile outf,EmbPMatMatch match, +// static AjBool silent_checktrans(AjPStr seq,AjPFile outf,EmbPMatMatch match, +static AjBool silent_checktrans(AjPStr seq,EmbPMatMatch match, AjPRinfo rlp, ajint begin, ajint radj, AjBool rev, ajint end, AjPSilent *res); -static void silent_fmt_sequence(AjPStr seq, AjPFile outf, ajint start, +/* static void silent_fmt_sequence(AjPStr seq, AjPFile outf, ajint start, AjBool num); -static void silent_fmt_hits(AjPList hits, AjPFile outf); +static void silent_fmt_hits(AjPList hits, AjPFile outf); */ static void silent_split_hits(AjPList *hits, AjPList *silents, AjPList *nonsilents, AjBool allmut); static ajint silent_basecompare(const void *a, const void *b); +static void silent_free(AjPSilent* sil); - - - +static void silent_make_report(AjPList hits, AjPReport report,AjPSeq seq, + AjBool compl); +static void silent_print_seq (AjPSeqout seqout,AjPStr seqname,AjPStr seq, + AjPStr extra); @@ -83,11 +87,13 @@ int main(int argc, char **argv) { AjPSeq seq=NULL; /*3 things we need from ACD/embInit*/ - AjPFile outf=NULL; + AjPSeqout outseq = NULL; /* sequence file */ + AjPReport report = NULL; /* enzyme report file */ AjPStr sstr=NULL; /*ajPatBruteForce requires the sequence as an AjPStr so we need one*/ AjPStr sname=NULL; /*and it wants a name as well*/ AjPStr revcomp=NULL; /*to hold reverse complement of sequence*/ + AjPStr tmpname=NULL; ajint RStotal; AjPStr enzymes=NULL; /* string for RE selection */ @@ -99,6 +105,7 @@ AjBool sshow; AjBool tshow; AjBool allmut; + AjBool comp; AjPList results1=NULL; /* for forward strand */ AjPList results2=NULL; /* for reverse strand */ @@ -111,10 +118,13 @@ /*and recover it from somewhere in memory*/ seq = ajAcdGetSeq("seq"); enzymes=ajAcdGetString("enzymes"); /* enzymes specified by user */ - sshow=ajAcdGetBool("sshow"); - tshow=ajAcdGetBool("tshow"); + report = ajAcdGetReport("outrep"); /* output report */ + outseq = ajAcdGetSeqout("outseq"); /* sequences output */ allmut=ajAcdGetBool("allmut"); - outf = ajAcdGetOutfile("outf"); + sshow=ajAcdGetBool("sshow"); + /* tshow = ajAcdGetBool("tshow"); /* display translated seq */ + /* Unable to make seqout to not make file if empty, so need sequence */ + tshow = ajTrue; shits=ajListNew(); nshits=ajListNew(); @@ -147,55 +157,59 @@ start=begin+1; if(sshow) { - ajFmtPrintF(outf,"SEQUENCE:\n"); - silent_fmt_sequence(sstr,outf,start,ajTrue); + silent_print_seq (outseq,sname,sstr,tmpname); } - results1=silent_mismatch(sstr,relist,outf,sname,RStotal,begin,radj, + results1=silent_mismatch(sstr,relist,outseq,sname,RStotal,begin,radj, ajFalse,end,tshow); silent_split_hits(&results1,&shits,&nshits,allmut); - ajFmtPrintF(outf,"Results for %S:\n\n",sname); - ajFmtPrintF(outf,"KEY:\n\tEnzyme\t\tEnzyme name\n" - "\tRS-Pattern\tRestriction enzyme recognition site pattern\n" - "\tMatch-Posn\tPosition of the first base of RS pattern in sequence\n" - "\tAA\t\tAmino acid. Original sequence(.)After mutation\n" - "\tBase-Posn\tPosition of base to be mutated in sequence\n" - "\tMutation\tThe base mutation to perform\n\n"); - - ajFmtPrintF(outf,"Silent mutations\n\n"); - silent_fmt_hits(shits,outf); - ajFmtPrintF(outf,"\n\n"); + ajReportSetHeader (report,tmpname); + silent_make_report (shits,report,seq,comp); + if(allmut) { - ajFmtPrintF(outf,"Non-Silent mutations\n\n"); + /* ajFmtPrintF(outf,"Non-Silent mutations\n\n"); silent_fmt_hits(nshits,outf); - ajFmtPrintF(outf,"\n\n"); + ajFmtPrintF(outf,"\n\n"); */ + ajStrAssC(&tmpname,"Non-Silent mutations\n"); + silent_make_report (nshits,report,seq,comp); } - if(sshow) + if (comp) { - ajFmtPrintF(outf,"REVERSE SEQUENCE:\n"); - silent_fmt_sequence(revcomp,outf,start,ajTrue); - } + if(sshow) + { + /* ajFmtPrintF(outf,"REVERSE SEQUENCE:\n"); + silent_fmt_sequence(revcomp,outf,start,ajTrue); */ + ajStrAssC (&tmpname, "Reverse complement sequence"); + silent_print_seq (outseq,sname,sstr,tmpname); + } - results2=silent_mismatch(revcomp,relist,outf,sname,RStotal,begin,radj, + results2=silent_mismatch(revcomp,relist,outseq,sname,RStotal,begin,radj, ajTrue,end,tshow); - silent_split_hits(&results2,&shits,&nshits,allmut); + silent_split_hits(&results2,&shits,&nshits,allmut); - ajFmtPrintF(outf,"\nResults for reverse of %S:\n\n",sname); - ajFmtPrintF(outf,"Silent mutations\n\n"); - silent_fmt_hits(shits,outf); - ajFmtPrintF(outf,"\n\n"); - if(allmut) - { - ajFmtPrintF(outf,"Non-Silent mutations\n\n"); - silent_fmt_hits(nshits,outf); - ajFmtPrintF(outf,"\n\n"); - } + /*ajFmtPrintF(outf,"\nResults for reverse of %S:\n\n",sname); + ajFmtPrintF(outf,"Silent mutations\n\n"); + silent_fmt_hits(shits,outf); + ajFmtPrintF(outf,"\n\n"); */ + ajStrAssC(&tmpname,"Reverse complement silent mutations\n"); + silent_make_report (shits,report,seq,comp); + if(allmut) + { + /* ajFmtPrintF(outf,"Non-Silent mutations\n\n"); + silent_fmt_hits(nshits,outf); + ajFmtPrintF(outf,"\n\n"); */ + ajStrAssC(&tmpname,"Revere complement non-silent mutations\n"); + silent_make_report (nshits,report,seq,comp); + } + } + ajReportClose (report); + ajReportDel (&report); ajSeqDel(&seq); ajStrDel(&revcomp); ajStrDel(&enzymes); @@ -205,7 +219,6 @@ ajListDel(&shits); ajListDel(&nshits); - ajFileClose(&outf); ajExit(); return 0; } @@ -235,7 +248,8 @@ ******************************************************************************/ -static AjPList silent_mismatch(AjPStr sstr, AjPList relist, AjPFile outf, +// static AjPList silent_mismatch(AjPStr sstr, AjPList relist, AjPFile outf, +static AjPList silent_mismatch(AjPStr sstr, AjPList relist, AjPSeqout outseq, AjPStr sname, ajint RStotal, ajint begin, ajint radj,AjBool rev, ajint end, AjBool tshow) @@ -244,6 +258,7 @@ AjPList results; AjPStr str; /*holds RS patterns*/ AjPStr tstr; + AjPStr tmp=NULL; AjBool dummy=ajFalse; /*need a bool for ajPatClassify*/ ajint mm; /*number of mismatches*/ ajint hits; /*number of pattern matches found*/ @@ -269,9 +284,11 @@ ajTrnStrFrame(table,sstr,1,&pep); /*1 stands for frame number*/ if(tshow) { - ajFmtPrintF(outf,"\n\nTRANSLATED SEQUENCE:\n"); + /* ajFmtPrintF(outf,"\n\nTRANSLATED SEQUENCE:\n"); silent_fmt_sequence(pep,outf,start,ajFalse); - ajFmtPrintF(outf,"\n\n"); + ajFmtPrintF(outf,"\n\n"); */ + ajStrAssC (&tmp," translated"); + silent_print_seq (outseq,sname,pep,tmp); } } else /* reverse strand */ @@ -280,9 +297,11 @@ ajTrnStrFrame(table,tstr,1,&pep); if(tshow) { - ajFmtPrintF(outf,"\n\nREVERSE TRANSLATED SEQUENCE:\n"); + /* ajFmtPrintF(outf,"\n\nREVERSE TRANSLATED SEQUENCE:\n"); silent_fmt_sequence(pep,outf,start,ajFalse); - ajFmtPrintF(outf,"\n\n"); + ajFmtPrintF(outf,"\n\n"); */ + ajStrAssC (&tmp," reverse complement translated"); + silent_print_seq (outseq,sname,pep,tmp); } } @@ -334,7 +353,7 @@ for(i=0;iissilent) { - ajListPushApp(*silents,(void *)res); - continue; + ajListPushApp(*silents,(void *)res); + continue; } if(allmut) - ajListPushApp(*nonsilents,(void *)res); + ajListPushApp(*nonsilents,(void *)res); else - { - ajStrDel(&res->code); - ajStrDel(&res->site); - ajStrDel(&res->seqaa); - ajStrDel(&res->reaa); - AJFREE(res); - } + silent_free (&res); } return; } +/* @funcstatic silent_make_report ******************************************** +** +** Print report +** +** @param [r] hits [AjPList] List of derestricted sites +** @param [r] report [AjPReport] Report to fill +** @param [r] seq [AjPSeq] Sequence +** @param [r] compl [AjBool] True if complementary strand +** @return [void] +******************************************************************************/ +static void silent_make_report(AjPList hits, AjPReport report,AjPSeq seq, + AjBool compl) +{ + AjPFeattable ftable = NULL; + AjPFeature seqf = NULL; + AjPStr tmp = NULL; + ajint len; + AjPSilent res; + + ajListSort(hits,silent_basecompare); + ftable = ajFeattableNewDna (ajSeqStr(seq)); + + while(ajListPop(hits,(void **)&res)) + { + len = ajStrLen (res->site); + if (compl) + seqf = ajFeatNew(ftable,NULL,NULL,res->match,res->match+len-1,0,'-',0); + else + seqf = ajFeatNewII (ftable,res->match,res->match+len-1); + ajFmtPrintS (&tmp,"*enzyme %S",res->code); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*site %S",res->site); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*origaa %S",res->seqaa); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*mutaa %S",res->reaa); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*base %d",res->base); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*orign %c",res->obase); + ajFeatTagAdd (seqf,NULL,tmp); + ajFmtPrintS (&tmp,"*mutn %c",res->nbase); + ajFeatTagAdd (seqf,NULL,tmp); + silent_free(&res); + } + + (void) ajReportWrite (report,ftable,seq); + ajFeatDel (&seqf); + ajFeattableDel (&ftable); + + return; +} + + +/* @funcstatic silent_print_seq ******************************************** +** +** Prints sequence +** +** @param [r] seqout [AjPSeqout] Sequence output +** @param [r] seqname [AjPStr] Sequence name +** @param [r] seq [AjPStr] Sequence +** @param [r] extra [AjPStr] Add to sequence name (can be NULL) +** @return [void] +******************************************************************************/ +static void silent_print_seq (AjPSeqout seqout,AjPStr seqname,AjPStr seq, + AjPStr extra) +{ + AjPSeq seqo = NULL; + AjPStr name = NULL; + seqo = ajSeqNew (); + ajStrAss (&name,seqname); + if (extra) + ajStrApp (&name,extra); + ajSeqAssName (seqo,name); + ajSeqAssSeq (seqo,seq); + ajSeqWrite (seqout,seqo); + + ajStrDel (&name); + ajSeqDel (&seqo); + + return; +} @@ -767,3 +863,21 @@ return((*(AjPSilent *)a)->base)-((*(AjPSilent *)b)->base); } + +/* @funcstatic silent_free ************************************************ +** +** Free allocated memory for silent mutant structure +** +** @param [d] sil [AjPSilent*] Silent mutant structure to be deleted +** @return [void] +******************************************************************************/ +static void silent_free(AjPSilent* sil) +{ + ajStrDel(&(*sil)->code); + ajStrDel(&(*sil)->site); + ajStrDel(&(*sil)->seqaa); + ajStrDel(&(*sil)->reaa); + AJFREE(sil); + + return; +}