1 /* @source banana application
2 **
3 ** Banana displays Bending and Curvature Calculations.
4 **
5 ** @author Copyright (C) Ian Longden (il@sanger.ac.uk)
6 ** @@
7 ** please reference the following report in any publication resulting from
8 ** use of this program.
9 ** Goodsell, D.S. & Dickerson, R.E. (1994) "Bending and Curvature Calculations
10 ** in B-DNA" Nucl. Acids. Res. 22, 5497-5503.
11 **
12 ** This program is free software; you can redistribute it and/or
13 ** modify it under the terms of the GNU General Public License
14 ** as published by the Free Software Foundation; either version 2
15 ** of the License, or (at your option) any later version.
16 **
17 ** This program is distributed in the hope that it will be useful,
18 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
19 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 ** GNU General Public License for more details.
21 **
22 ** You should have received a copy of the GNU General Public License
23 ** along with this program; if not, write to the Free Software
24 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
25 ******************************************************************************/
26 #include "emboss.h"
27 #include <math.h>
28
29
30 #define sinfban(x) (float)sin((double)x)
31 #define cosfban(x) (float)cos((double)x)
32
33
34
35
36 /* @prog banana ***************************************************************
37 **
38 ** Bending and curvature plot in B-DNA
39 **
40 ******************************************************************************/
41
main(int argc,char ** argv)42 int main(int argc, char **argv)
43 {
44 AjPSeq seq;
45 AjPGraph graph = 0;
46 AjPFile outf = NULL;
47 AjPFile file = NULL;
48 AjPStr buffer = NULL;
49 float twist[4][4][4];
50 float roll[4][4][4];
51 float tilt[4][4][4];
52 float rbend;
53 float rcurve;
54 float bendscale;
55 float curvescale;
56 float twistsum = (float) 0.0;
57 float pi = (float) 3.14159;
58 float pifac = (pi/(float) 180.0);
59 float pi2 = pi/(float) 2.0;
60 ajint *iseq = NULL;
61 float *x;
62 float *y;
63 float *xave;
64 float *yave;
65 float *curve;
66 float *bend;
67 const char *ptr;
68 ajint i;
69 ajint ii;
70 ajint k;
71 ajint j;
72 char residue[2];
73 float maxbend;
74 float minbend;
75 float bendfactor;
76 float maxcurve;
77 float mincurve;
78 float curvefactor;
79
80 float fxp;
81 float fyp;
82 float yincr;
83 float yy1;
84 ajint ixlen;
85 ajint iylen;
86 ajint ixoff;
87 ajint iyoff;
88 float ystart;
89 float defheight = 0.0;
90 float currentheight = 0.0;
91 ajint count;
92 ajint portrait = 0;
93 ajint title = 0;
94 ajint numres;
95 ajint ibeg;
96 ajint iend;
97 ajint ilen;
98 AjPStr sstr = NULL;
99 ajint ipos;
100 float dx;
101 float dy;
102 float rxsum;
103 float rysum;
104 float yp1;
105 float yp2;
106 double td;
107
108 embInit("banana", argc, argv);
109 seq = ajAcdGetSeq("sequence");
110 file = ajAcdGetDatafile("anglesfile");
111 outf = ajAcdGetOutfile("outfile");
112 graph = ajAcdGetGraph("graph");
113 numres = ajAcdGetInt("residuesperline");
114
115 ibeg = ajSeqGetBegin(seq);
116 iend = ajSeqGetEnd(seq);
117
118 ajStrAssignSubS(&sstr, ajSeqGetSeqS(seq), ibeg-1, iend-1);
119 ilen = ajStrGetLen(sstr);
120
121 AJCNEW0(iseq,ilen+1);
122 AJCNEW0(x,ilen+1);
123 AJCNEW0(y,ilen+1);
124 AJCNEW0(xave,ilen+1);
125 AJCNEW0(yave,ilen+1);
126 AJCNEW0(curve,ilen+1);
127 AJCNEW0(bend,ilen+1);
128
129 ptr= ajStrGetPtr(sstr);
130
131 for(ii=0;ii<ilen;ii++)
132 {
133 if(*ptr=='A' || *ptr=='a')
134 iseq[ii+1] = 0;
135 else if(*ptr=='T' || *ptr=='t')
136 iseq[ii+1] = 1;
137 else if(*ptr=='G' || *ptr=='g')
138 iseq[ii+1] = 2;
139 else if(*ptr=='C' || *ptr=='c')
140 iseq[ii+1] = 3;
141 else
142 ajErr("%c is not an ATCG hence not valid",*ptr);
143
144 ptr++;
145 }
146
147
148 if(!file)
149 ajErr("Banana failed to open angle file");
150
151 ajReadline(file,&buffer); /* 3 junk lines */
152 ajReadline(file,&buffer);
153 ajReadline(file,&buffer);
154
155 for(k=0;k<4;k++)
156 for(ii=0;ii<4;ii++)
157 {
158 if(ajReadline(file,&buffer))
159 {
160 sscanf(ajStrGetPtr(buffer),"%f,%f,%f,%f",
161 &twist[ii][0][k],&twist[ii][1][k],&twist[ii][2][k],
162 &twist[ii][3][k]);
163 }
164 else
165 ajErr("Error reading angle file");
166
167 for(j=0;j<4;j++)
168 twist[ii][j][k] *= pifac;
169 }
170
171
172 for(k=0;k<4;k++)
173 for(ii=0;ii<4;ii++)
174 if(ajReadline(file,&buffer))
175 {
176 sscanf(ajStrGetPtr(buffer),"%f,%f,%f,%f",&roll[ii][0][k],
177 &roll[ii][1][k],&roll[ii][2][k],&roll[ii][3][k]);
178 }
179 else
180 ajErr("Error reading angle file");
181
182
183 for(k=0;k<4;k++)
184 for(ii=0;ii<4;ii++)
185 if(ajReadline(file,&buffer))
186 sscanf(ajStrGetPtr(buffer),"%f,%f,%f,%f",&tilt[ii][0][k],
187 &tilt[ii][1][k],&tilt[ii][2][k],&tilt[ii][3][k]);
188 else
189 ajErr("Error reading angle file");
190
191
192 if(ajReadline(file,&buffer))
193 sscanf(ajStrGetPtr(buffer),"%f,%f,%f,%f",&rbend,&rcurve,
194 &bendscale,&curvescale);
195 else
196 ajErr("Error reading angle file");
197
198 ajFileClose(&file);
199 ajStrDel(&buffer);
200
201
202 for(ii=1;ii<ilen-1;ii++)
203 {
204 twistsum += twist[iseq[ii]][iseq[ii+1]][iseq[ii+2]];
205 dx = (roll[iseq[ii]][iseq[ii+1]][iseq[ii+2]]*sinfban(twistsum)) +
206 (tilt[iseq[ii]][iseq[ii+1]][iseq[ii+2]]*sinfban(twistsum-pi2));
207 dy = roll[iseq[ii]][iseq[ii+1]][iseq[ii+2]]*cosfban(twistsum) +
208 tilt[iseq[ii]][iseq[ii+1]][iseq[ii+2]]*cosfban(twistsum-pi2);
209
210 x[ii+1] = x[ii]+dx;
211 y[ii+1] = y[ii]+dy;
212
213 }
214
215 for(ii=6;ii<ilen-6;ii++)
216 {
217 rxsum = 0.0;
218 rysum = 0.0;
219 for(k=-4;k<=4;k++)
220 {
221 rxsum+=x[ii+k];
222 rysum+=y[ii+k];
223 }
224 rxsum+=(x[ii+5]*(float)0.5);
225 rysum+=(y[ii+5]*(float)0.5);
226 rxsum+=(x[ii-5]*(float)0.5);
227 rysum+=(y[ii-5]*(float)0.5);
228
229 xave[ii] = rxsum*(float)0.1;
230 yave[ii] = rysum*(float)0.1;
231 }
232
233 for(i=(ajint)rbend+1;i<=ilen-(ajint)rbend-1;i++)
234 {
235 td = sqrt(((x[i+(ajint)rbend]-x[i-(ajint)rbend])*
236 (x[i+(ajint)rbend]-x[i-(ajint)rbend])) +
237 ((y[i+(ajint)rbend]-y[i-(ajint)rbend])*
238 (y[i+(ajint)rbend]-y[i-(ajint)rbend])));
239 bend[i] = (float) td;
240 bend[i]*=bendscale;
241 }
242
243 for(i=(ajint)rcurve+6;i<=ilen-(ajint)rcurve-6;i++)
244 {
245 td = sqrt(((xave[i+(ajint)rcurve]-
246 xave[i-(ajint)rcurve])*(xave[i+(ajint)rcurve]-
247 xave[i-(ajint)rcurve]))+
248 ((yave[i+(ajint)rcurve]-yave[i-(ajint)rcurve])*
249 (yave[i+(ajint)rcurve]-yave[i-(ajint)rcurve])));
250 curve[i] = (float) td;
251 }
252
253
254 if(outf)
255 {
256 ajFmtPrintF(outf,"Base Bend Curve\n");
257 ptr = ajStrGetPtr(sstr);
258 for(ii=1;ii<=ilen;ii++)
259 {
260 ajFmtPrintF(outf,"%c %6.1f %6.1f\n",
261 *ptr, bend[ii], curve[ii]);
262 ptr++;
263 }
264 ajFileClose(&outf);
265 }
266
267 if(graph)
268 {
269 maxbend = minbend = 0.0;
270 maxcurve = mincurve = 0.0;
271 for(ii=1;ii<=ilen;ii++)
272 {
273 if(bend[ii] > maxbend)
274 maxbend = bend[ii];
275 if(bend[ii] < minbend)
276 minbend = bend[ii];
277 if(curve[ii] > maxcurve)
278 maxcurve = curve[ii];
279 if(curve[ii] < mincurve)
280 mincurve = curve[ii];
281 }
282
283 ystart = 75.0;
284
285 ajGraphAppendTitleS(graph, ajSeqGetUsaS(seq));
286
287 ajGraphicsSetPagesize(960, 768);
288
289 ajGraphOpenWin(graph,(float)-1.0, (float)numres+(float)10.0,
290 (float)0.0, ystart+(float)5.0);
291
292 ajGraphicsGetParamsPage(&fxp,&fyp,&ixlen,&iylen,&ixoff,&iyoff);
293
294 if(portrait)
295 {
296 ixlen = 768;
297 iylen = 960;
298 }
299 else
300 {
301 ixlen = 960;
302 iylen = 768;
303 }
304
305 ajGraphicsGetCharsize(&defheight,¤theight);
306 if(!currentheight)
307 {
308 defheight = currentheight = (float) 4.440072;
309 currentheight = defheight *
310 ((float)ixlen/ ((float)(numres)*(currentheight+(float)1.0)))
311 /currentheight;
312 }
313 ajGraphicsSetCharscale(((float)ixlen/((float)(numres)*
314 (currentheight+(float)1.0)))/
315 currentheight);
316 ajGraphicsGetCharsize(&defheight,¤theight);
317
318 yincr = (currentheight + (float)3.0)*(float)0.3;
319
320 if(!title)
321 yy1 = ystart;
322 else
323 yy1 = ystart-(float)5.0;
324
325 count = 1;
326
327 residue[1]='\0';
328
329 bendfactor = (3*yincr)/maxbend;
330 curvefactor = (3*yincr)/maxcurve;
331
332 ptr = ajStrGetPtr(sstr);
333
334 yy1 = yy1-(yincr*((float)5.0));
335 for(ii=1;ii<=ilen;ii++)
336 {
337 if(count > numres)
338 {
339 yy1 = yy1-(yincr*((float)5.0));
340 if(yy1<1.0)
341 {
342 if(!title)
343 yy1=ystart;
344 else
345 yy1 = ystart-(float)5.0;
346
347 yy1 = yy1-(yincr*((float)5.0));
348 ajGraphNewpage(graph,AJFALSE);
349 }
350 count = 1;
351 }
352 residue[0] = *ptr;
353
354 ajGraphicsDrawposTextAtend((float)(count)+(float)2.0,yy1,residue);
355
356 if(ii>1 && ii < ilen)
357 {
358 yp1 = yy1+yincr + (bend[ii]*bendfactor);
359 yp2 = yy1+yincr + (bend[ii+1]*bendfactor);
360 ajGraphicsDrawposLine((float)count+(float)1.5,yp1,
361 (float)(count)+(float)2.5,yp2);
362 }
363
364 ipos = ilen-(ajint)rcurve-7;
365 if(ipos < 0)
366 ipos = 0;
367
368 if(ii>rcurve+5 && ii<ipos)
369 {
370 yp1 = yy1+yincr + (curve[ii]*curvefactor);
371 yp2 = yy1+yincr + (curve[ii+1]*curvefactor);
372 ajGraphicsDrawposLine((float)count+(float)1.7,yp1,
373 (float)(count)+(float)2.3,yp2);
374 }
375
376 ajGraphicsDrawposLine((float)count+(float)1.5,yy1+yincr,
377 (float)(count)+(float)2.5,yy1+yincr);
378
379 count++;
380 ptr++;
381 }
382
383 ajGraphicsClose();
384 }
385
386 AJFREE(iseq);
387 AJFREE(x);
388 AJFREE(y);
389 AJFREE(xave);
390 AJFREE(yave);
391 AJFREE(curve);
392 AJFREE(bend);
393
394 ajStrDel(&sstr);
395
396 ajSeqDel(&seq);
397 ajFileClose(&file);
398 ajFileClose(&outf);
399 ajGraphxyDel(&graph);
400
401 embExit();
402
403 return 0;
404 }
405