/* @source profit application ** ** Scan a protein database or sequence with a profile or matrix ** @author Copyright (C) Alan Bleasby (ableasby@hgmp.mrc.ac.uk) ** @@ ** ** This program is free software; you can redistribute it and/or ** modify it under the terms of the GNU General Public License ** as published by the Free Software Foundation; either version 2 ** of the License, or (at your option) any later version. ** ** This program is distributed in the hope that it will be useful, ** but WITHOUT ANY WARRANTY; without even the implied warranty of ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ** GNU General Public License for more details. ** ** You should have received a copy of the GNU General Public License ** along with this program; if not, write to the Free Software ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ******************************************************************************/ #include "emboss.h" #include #define AZ 27 float maxscore = 0.0; ajuint seqnum=0; static void profit_read_profile(AjPFile inf, AjPStr *name, AjPStr *mname, ajint *mlen, float *gapopen, float *gapextend, ajint *thresh, float *maxs, AjPStr *cons); static void profit_read_simple(AjPFile inf, AjPStr *name, ajint *mlen, ajint *maxs, ajint *thresh,AjPStr *cons); static void profit_scan_profile(const AjPStr substr, const AjPStr seqname, ajint mlen, float * const *fmatrix, ajint thresh, float maxs, AjPFile outf); static void profit_scan_simple(const AjPStr substr, const AjPStr seqname, ajint mlen, ajint maxs, ajint thresh, ajint * const *matrix, AjPFile outf); static void profit_printHits(const AjPStr seqname, ajint pos, ajint score, AjPFile outf); static ajint profit_getType(AjPFile inf); /* @prog profit *************************************************************** ** ** Scan a sequence or database with a matrix or profile ** ******************************************************************************/ int main(int argc, char **argv) { AjPSeqall seqall; AjPSeq seq = NULL; AjPFile inf = NULL; AjPFile outf = NULL; AjPStr strand = NULL; AjPStr substr = NULL; AjPStr name = NULL; AjPStr mname = NULL; AjPStr seqname = NULL; AjPStr line = NULL; AjPStr cons = NULL; ajint type; ajint begin; ajint end; ajint i; ajint j; ajint **matrix = NULL; float **fmatrix = NULL; void **fptr = NULL; ajint mlen; ajint maxs; float maxfs; ajint thresh; float gapopen; float gapextend; const char *p; embInit("profit", argc, argv); seqall = ajAcdGetSeqall("sequence"); inf = ajAcdGetInfile("infile"); outf = ajAcdGetOutfile("outfile"); substr = ajStrNew(); name = ajStrNew(); mname = ajStrNew(); line = ajStrNew(); if(!(type=profit_getType(inf))) ajFatal("Unrecognised profile/matrix file format"); switch(type) { case 1: profit_read_simple(inf, &name, &mlen, &maxs, &thresh, &cons); AJCNEW(matrix, mlen); fptr=(void **)matrix; for(i=0;i= threshold %d (max score %d)\n#\n",thresh, maxs); break; case 2: profit_read_profile(inf,&name,&mname,&mlen,&gapopen,&gapextend,&thresh, &maxfs, &cons); AJCNEW(fmatrix, mlen); fptr=(void **)fmatrix; for(i=0;i= threshold %d (max score %.2f)\n#\n", thresh,maxfs); break; case 3: profit_read_profile(inf,&name,&mname,&mlen,&gapopen,&gapextend,&thresh, &maxfs, &cons); AJCNEW(fmatrix, mlen); fptr=(void **)fmatrix; for(i=0;i= threshold %d (max score %.2f)\n#\n", thresh,maxfs); break; default: ajFatal("Switch type error"); break; } while(ajSeqallNext(seqall, &seq)) { seqnum++; begin = ajSeqallGetseqBegin(seqall); end = ajSeqallGetseqEnd(seqall); ajStrAssignC(&seqname,ajSeqGetNameC(seq)); strand = ajSeqGetSeqCopyS(seq); ajStrAssignSubC(&substr,ajStrGetPtr(strand),begin-1,end-1); switch(type) { case 1: profit_scan_simple(substr,seqname,mlen,maxs, thresh,matrix,outf); break; case 2: profit_scan_profile(substr,seqname,mlen,fmatrix, thresh,maxfs,outf); break; case 3: profit_scan_profile(substr,seqname,mlen,fmatrix, thresh,maxfs,outf); break; default: break; } ajStrDel(&strand); } for(i=0;i= thresh) profit_printHits(seqname,i,score,outf); } return; } /* @funcstatic profit_printHits *********************************************** ** ** Print results for profit ** ** @param [r] seqname [const AjPStr] sequence name ** @param [r] pos [ajint] position ** @param [r] score [ajint] score ** @param [u] outf [AjPFile] outfile ** @@ ******************************************************************************/ static void profit_printHits(const AjPStr seqname, ajint pos, ajint score, AjPFile outf) { ajFmtPrintF(outf,"%S %d Percentage: %d\n",seqname,pos+1,score); return; } /* @funcstatic profit_scan_profile ******************************************** ** ** Scan sequence with a profile ** ** @param [r] substr [const AjPStr] sequence ** @param [r] seqname [const AjPStr] sequence name ** @param [r] mlen [ajint] profile length ** @param [r] fmatrix [float* const *] score matrix ** @param [r] thresh [ajint] threshold ** @param [r] maxs [float] max score ** @param [u] outf [AjPFile] outfile ** @@ ******************************************************************************/ static void profit_scan_profile (const AjPStr substr, const AjPStr seqname, ajint mlen, float * const *fmatrix, ajint thresh, float maxs, AjPFile outf) { ajint len; ajint i; ajint j; ajint lim; const char *p; float score; float sum; len = ajStrGetLen(substr); lim = len-mlen+1; p = ajStrGetPtr(substr); ajUser("scan_profile mlen:%d len:%d lim:%d '%S'", mlen, len, lim, seqname); for(i=0;i= thresh) profit_printHits(seqname,i,(ajint)score,outf); if(score > maxscore) { ajUser("maxscore: %f sum: %f i:%d '%S' %u maxs:%f", score, sum, i, seqname, seqnum, maxs); maxscore = score; } } return; }