1 /* @source hmoment application
2 **
3 ** Calculate hydrophobic moment
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
26
27
28
29 static float hmoment_calchm(const char *p, int pos, int window, ajint angle);
30 static void hmoment_addgraph(AjPGraph graph, ajint limit,
31 const float *x, const float *y,
32 float ymax, ajint colour, ajint angle,
33 ajint window, float baseline, const char *sname);
34
35
36
37
38 /* @prog hmoment **************************************************************
39 **
40 ** Hydrophobic moment calculation
41 **
42 ******************************************************************************/
43
main(int argc,char ** argv)44 int main(int argc, char **argv)
45 {
46 AjPSeqall seqall;
47 AjPSeq seq;
48 AjPFile outf;
49 AjPStr str = NULL;
50 AjPStr st = NULL;
51
52 AjBool plot;
53 AjPGraph graph = NULL;
54 AjBool twin;
55 float baseline;
56 float ymax;
57
58 ajint beg;
59 ajint end;
60 ajint window;
61 ajint len;
62 ajint limit;
63
64 ajint aangle;
65 ajint bangle;
66
67 ajint i;
68
69 float *x = NULL;
70 float *ya = NULL;
71 float *yb = NULL;
72
73 const char *p;
74 const char *sname;
75
76
77 embInit("hmoment", argc, argv);
78
79 seqall = ajAcdGetSeqall("seqall");
80 plot = ajAcdGetToggle("plot");
81 window = ajAcdGetInt("window");
82 aangle = ajAcdGetInt("aangle");
83 bangle = ajAcdGetInt("bangle");
84 baseline = ajAcdGetFloat("baseline");
85 twin = ajAcdGetBoolean("double");
86
87 /* only one will be used - see variable 'plot' */
88
89 outf = ajAcdGetOutfile("outfile");
90 graph = ajAcdGetGraphxy("graph");
91
92
93 str = ajStrNew();
94 st = ajStrNew();
95
96 while(ajSeqallNext(seqall, &seq))
97 {
98 beg = ajSeqallGetseqBegin(seqall);
99 end = ajSeqallGetseqEnd(seqall);
100 len = end-beg+1;
101
102 limit = len-window+1;
103 sname = ajSeqGetNameC(seq);
104
105 ajStrAssignSubC(&str,ajSeqGetSeqC(seq),--beg,--end);
106 ajStrFmtUpper(&str);
107 p = ajStrGetPtr(str);
108
109 if(limit>0)
110 {
111 AJCNEW0(x,limit);
112 AJCNEW0(ya,limit);
113 AJCNEW0(yb,limit);
114 }
115
116 for(i=0;i<limit;++i)
117 {
118 x[i] = (float)i+1+beg;
119 ya[i] = hmoment_calchm(p,i,window,aangle);
120 yb[i] = hmoment_calchm(p,i,window,bangle);
121 }
122
123 for(i=0,ymax=-100.;i<limit;++i)
124 {
125 ymax = ymax > ya[i] ? ymax : ya[i];
126 if(twin)
127 ymax = ymax > yb[i] ? ymax : yb[i];
128 }
129
130 if(!plot && outf)
131 {
132 ajFmtPrintF(outf,"HMOMENT of %s from %d to %d\n\n",sname,
133 beg+1,end+1);
134 if(twin)
135 {
136 ajFmtPrintF(outf,"Window: %d First angle: %d second angle:"
137 " %d Max uH: %.3f\n",window, aangle, bangle,
138 ymax);
139 ajFmtPrintF(outf,"Position\tFirst\tSecond\n");
140 for(i=0;i<limit;++i)
141 ajFmtPrintF(outf,"%d\t\t%.3lf\t%.3lf\n",(int)x[i],ya[i],
142 yb[i]);
143 }
144 else
145 {
146 ajFmtPrintF(outf,"Window: %d Angle: %d Max uH: %.3f\n",
147 window, aangle, ymax);
148 ajFmtPrintF(outf,"Position\tuH\n");
149 for(i=0;i<limit;++i)
150 ajFmtPrintF(outf,"%d\t\t%.3lf\n",(int)x[i],ya[i]);
151 }
152
153 }
154 else if (plot)
155 {
156 if(twin)
157 ajGraphSetMulti(graph,2);
158 else
159 ajGraphSetMulti(graph,1);
160
161 ajGraphxySetflagOverlay(graph,ajFalse);
162 ajGraphSetDatanameS(graph, ajSeqGetNameS(seq));
163
164 hmoment_addgraph(graph,limit,x,ya,ymax,BLACK,aangle,window,
165 baseline, sname);
166 if(twin)
167 hmoment_addgraph(graph,limit,x,yb,ymax,RED,bangle,window,
168 baseline, sname);
169
170 if(limit>0)
171 ajGraphxyDisplay(graph,ajFalse);
172 }
173
174 if(limit>0)
175 {
176 AJFREE(x);
177 AJFREE(ya);
178 AJFREE(yb);
179 }
180 }
181
182 ajGraphicsClose();
183 ajGraphxyDel(&graph);
184 ajFileClose(&outf);
185
186 ajStrDel(&str);
187 ajStrDel(&st);
188
189 ajSeqallDel(&seqall);
190 ajSeqDel(&seq);
191 embExit();
192
193 return 0;
194 }
195
196
197
198
199 /* @funcstatic hmoment_addgraph ***********************************************
200 **
201 ** Undocumented.
202 **
203 ** @param [u] graph [AjPGraph] Undocumented
204 ** @param [r] limit [ajint] Undocumented
205 ** @param [r] x [const float*] Undocumented
206 ** @param [r] y [const float*] Undocumented
207 ** @param [r] ymax [float] Undocumented
208 ** @param [r] colour [ajint] Undocumented
209 ** @param [r] angle [ajint] Undocumented
210 ** @param [r] window [ajint] Undocumented
211 ** @param [r] baseline [float] Undocumented
212 ** @param [r] sname [const char*] Sequence name
213 ** @@
214 ******************************************************************************/
215
216
hmoment_addgraph(AjPGraph graph,ajint limit,const float * x,const float * y,float ymax,ajint colour,ajint angle,ajint window,float baseline,const char * sname)217 static void hmoment_addgraph(AjPGraph graph, ajint limit,
218 const float *x, const float *y,
219 float ymax, ajint colour, ajint angle,
220 ajint window, float baseline, const char *sname)
221 {
222 ajint i;
223
224 AjPGraphdata data;
225 AjPStr st = NULL;
226
227 if(limit<1)
228 return;
229
230
231 data = ajGraphdataNewI(limit);
232
233 st = ajStrNew();
234
235 for(i=0;i<limit;++i)
236 {
237 data->x[i] = x[i];
238 data->y[i] = y[i];
239 }
240
241 ajGraphdataSetColour(data,colour);
242 ajGraphdataSetMinmax(data,x[0],x[limit-1],0.,ymax);
243 ajGraphdataSetTruescale(data,x[0],x[limit-1],0.,ymax);
244
245 ajGraphdataSetTypeC(data,"2D Plot Float");
246
247 ajFmtPrintS(&st,"HMOMENT of %s. Window:%d",sname,window);
248 ajGraphdataSetTitleS(data,st);
249
250 ajFmtPrintS(&st,"uH (%d deg)",angle);
251 ajGraphdataSetYlabelS(data,st);
252
253 ajFmtPrintS(&st,"Position (w=%d)",window);
254 ajGraphdataSetXlabelS(data,st);
255
256 ajGraphdataAddposLine(data,x[0],baseline,x[limit-1],baseline,BLUE);
257
258 ajGraphDataAdd(graph,data);
259
260 ajStrDel(&st);
261
262 return;
263 }
264
265
266
267
268 /* @funcstatic hmoment_calchm *************************************************
269 **
270 ** Undocumented.
271 **
272 ** @param [r] p [const char*] Undocumented
273 ** @param [r] pos [int] Undocumented
274 ** @param [r] window [int] Undocumented
275 ** @param [r] angle [ajint] Undocumented
276 ** @return [float] Undocumented
277 ** @@
278 ******************************************************************************/
279
280
hmoment_calchm(const char * p,int pos,int window,ajint angle)281 static float hmoment_calchm(const char *p, int pos, int window, ajint angle)
282 {
283 ajint i;
284 double h;
285 ajint res;
286 double sumsin;
287 double sumcos;
288 double tangle;
289 double hm;
290
291 /*
292 ** B interpolated from D and N using Dayhoff frequencies
293 ** Z interpolated from E and Q using Dayhoff frequencies
294 ** J interpolated from I and L using Dayhoff frequencies
295 ** X average of all using Dayhoff frequencies
296 ** O and U set to X due to lack of available data
297 */
298 double hydata[]=
299 {
300 0.62, -0.84, 0.29, -0.90, -0.74, 1.19, 0.48, -0.40, /* ABCDEFGH */
301 1.38, 1.18, -1.50, 1.06, 0.64, -0.78, -0.68, 0.12, /* IJKLMNOP */
302 -0.85, -2.53, -0.18, -0.05, -0.68, 1.08, 0.81, -0.68, /* QRSTUVWX */
303 0.26, -0.78 /* YZ */
304 };
305
306 sumsin = sumcos = (double)0.;
307 tangle = angle;
308
309 for(i=0;i<window;++i)
310 {
311 res = p[pos+i];
312 h = hydata[ajBasecodeToInt(res)];
313
314 sumsin += (h * sin(ajCvtDegToRad((float)tangle)));
315 sumcos += (h * cos(ajCvtDegToRad((float)tangle)));
316 tangle = (double) (((ajint)tangle+angle) % 360);
317 }
318
319
320 sumsin *= sumsin;
321 sumcos *= sumcos;
322
323
324 hm = sqrt(sumsin+sumcos) / (double)window;
325
326 return (float)hm;
327 }
328