1 /* @source freak application
2 **
3 ** Calculate residue frequencies using sliding window
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
25
26
27
28 /* @prog freak ****************************************************************
29 **
30 ** Residue/base frequency table or plot
31 **
32 ******************************************************************************/
33
main(int argc,char ** argv)34 int main(int argc, char **argv)
35 {
36 AjPSeqall seqall;
37 AjPSeq seq;
38 AjPFile outf;
39 AjPStr bases = NULL;
40 AjPStr str = NULL;
41 AjBool plot;
42 AjPGraph graph=NULL;
43 AjPGraphdata fgraph=NULL;
44 AjPStr st = NULL;
45
46 ajint c;
47 ajint pos;
48 ajint end;
49 ajint step;
50 ajint window;
51 ajint t;
52
53 ajint i;
54 ajint j;
55 ajint k;
56 const char *p;
57 const char *q;
58 float f;
59
60 float *x = NULL;
61 float *y = NULL;
62
63 float max = 0.;
64 float min = 0.;
65
66
67 embInit("freak", argc, argv);
68
69 seqall = ajAcdGetSeqall("seqall");
70 plot = ajAcdGetToggle("plot");
71 step = ajAcdGetInt("step");
72 window = ajAcdGetInt("window");
73 bases = ajAcdGetString("letters");
74
75 /* only one will be used - see variable 'plot' */
76
77 outf = ajAcdGetOutfile("outfile");
78 graph = ajAcdGetGraphxy("graph");
79
80
81 st = ajStrNew();
82
83 ajStrFmtUpper(&bases);
84 q = ajStrGetPtr(bases);
85
86 while(ajSeqallNext(seqall, &seq))
87 {
88 pos = ajSeqallGetseqBegin(seqall);
89 end = ajSeqallGetseqEnd(seqall);
90
91 str = ajSeqGetSeqCopyS(seq);
92 ajStrFmtUpper(&str);
93 p = ajStrGetPtr(str);
94
95 c = 0;
96 --pos;
97 --end;
98 t = pos;
99 while(t+window <= end+1)
100 {
101 ++c;
102 t += step;
103 }
104
105
106 if(c)
107 {
108 x = (float *) AJALLOC(c * sizeof(float));
109 y = (float *) AJALLOC(c * sizeof(float));
110 }
111
112
113 for(i=0;i<c;++i)
114 {
115 t = i*step+pos;
116 x[i] = (float)(t+1) /*+ window/2.0*/;
117 f = 0.;
118 for(j=0;j<window;++j)
119 {
120 k = t + j;
121 if(strchr(q,p[k]))
122 ++f;
123 }
124 y[i] = f / (float)window;
125 }
126
127 if(!plot && c && outf)
128 {
129 ajFmtPrintF(outf,"FREAK of %s from %d to %d Window %d Step %d\n\n",
130 ajSeqGetNameC(seq),pos+1,end+1,window,step);
131 for(i=0;i<c;++i)
132 ajFmtPrintF(outf,"%-10d %f\n",(ajint)x[i],y[i]);
133 }
134 else if(plot && c)
135 {
136 fgraph = ajGraphdataNewI(c);
137 ajGraphSetDatanameS(graph, ajSeqGetNameS(seq));
138 ajGraphSetTitleS(graph,ajSeqGetNameS(seq));
139 ajFmtPrintS(&st,"From %d to %d. Residues:%s Window:%d Step:%d",
140 pos+1,end+1,ajStrGetPtr(bases),window,step);
141 ajGraphSetSubtitleS(graph,st);
142 ajGraphSetXlabelC(graph,"Position");
143 ajGraphSetYlabelC(graph,"Frequency");
144 ajGraphxySetXstartF(graph,x[0]);
145 ajGraphxySetXendF(graph,x[c-1]);
146 ajGraphxySetYstartF(graph,0.);
147 ajGraphxySetYendF(graph,y[c-1]);
148 ajGraphxySetXrangeII(graph,(ajint)x[0],(ajint)x[c-1]);
149 ajGraphxySetYrangeII(graph,0,(ajint)y[c-1]);
150 ajGraphdataSetMinmax(fgraph,x[0],x[c-1],0.,1.0);
151 ajGraphicsCalcRange(y,c,&min,&max);
152 ajGraphdataSetTruescale(fgraph,x[0],x[c-1],min,max);
153 ajGraphdataSetTypeC(fgraph,"2D Plot");
154
155 ajGraphdataAddXY(fgraph,x,y);
156 ajGraphDataReplace(graph,fgraph);
157 ajGraphxyDisplay(graph,ajFalse);
158 }
159
160
161 AJFREE(x);
162 AJFREE(y);
163 ajStrDel(&str);
164 }
165
166 ajGraphicsClose();
167 ajGraphxyDel(&graph);
168
169 ajFileClose(&outf);
170 ajSeqallDel(&seqall);
171 ajSeqDel(&seq);
172
173 ajStrDel(&str);
174 ajStrDel(&bases);
175 ajStrDel(&st);
176
177 embExit();
178
179 return 0;
180 }
181