1 /* @source pepinfo application
2 **
3 ** Displays properties of the amino acid residues in a peptide sequence
4 ** @author Copyright (C) Mark Faller (mfaller@hgmp.mrc.ac.uk)
5 ** @@
6 **
7 ** This program is free software; you can redistribute it and/or
8 ** modify it under the terms of the GNU General Public License
9 ** as published by the Free Software Foundation; either version 2
10 ** of the License, or (at your option) any later version.
11 **
12 ** This program is distributed in the hope that it will be useful,
13 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 ** GNU General Public License for more details.
16 **
17 ** You should have received a copy of the GNU General Public License
18 ** along with this program; if not, write to the Free Software
19 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 ******************************************************************************/
21
22 #include "emboss.h"
23
24
25 #define NOY (AJGRAPH_X_BOTTOM + AJGRAPH_Y_LEFT + AJGRAPH_Y_RIGHT + \
26 AJGRAPH_Y_INVERT_TICK + AJGRAPH_X_INVERT_TICK + AJGRAPH_X_TICK + \
27 AJGRAPH_X_LABEL + AJGRAPH_Y_LABEL + AJGRAPH_TITLE)
28
29
30
31
32 static void pepinfo_plotHistInt2(AjPHist hist,
33 const ajint * results, ajint hist_num,
34 const char* header,
35 const char* xtext, const char * ytext);
36 static void pepinfo_plotGraph2Float(AjPGraph graphs,
37 const float * results,
38 const char* title_text,
39 const char * xtext, const char * ytext);
40 static void pepinfo_printFloatResults(AjPFile outfile, const AjPSeq seq,
41 const float * results,
42 const char* header);
43 static void pepinfo_printIntResults(AjPFile outfile, const AjPSeq seq,
44 const ajint* results, const char* header);
45
46
47
48
49 static ajint seq_start;
50 static ajint seq_end;
51 static ajint win_mid;
52 static ajint seq_begin;
53
54
55
56
57 /* @prog pepinfo **************************************************************
58 **
59 ** Plots simple amino acid properties in parallel
60 **
61 ******************************************************************************/
62
main(int argc,char ** argv)63 int main(int argc, char **argv)
64 {
65
66 AjPSeq inseq = NULL;
67 AjPFile outfile = NULL;
68 ajint hwindow;
69 AjPFile aa_properties = NULL;
70 AjPFile aa_hydropathy = NULL;
71
72 AjBool do_seq;
73 AjBool do_general;
74 AjBool do_hydropathy;
75 AjBool do_plot;
76 AjPStr key;
77 const AjPStr value;
78 AjIList listIter;
79 AjPTable table;
80
81 AjPStr tmpa = NULL;
82 AjPStr tmpb = NULL;
83 ajint i;
84 ajint j;
85 ajint k;
86 ajint cnt;
87
88 ajint * ival;
89 ajint * iv[9];
90 float * pfloat;
91 float * pf[3];
92 float num;
93 float total;
94
95 AjPGraph graphs = NULL;
96 AjPHist hist = NULL;
97 ajint numGraphs = 0;
98
99 /* Data_table aa_props, aa_hydro, aa_acc;*/
100 AjPList aa_props = NULL;
101 AjPList aa_hydro = NULL;
102
103 const char * propertyTitles[] =
104 {
105 "Tiny",
106 "Small",
107 "Aliphatic",
108 "Aromatic",
109 "Non-polar",
110 "Polar",
111 "Charged",
112 "Positive",
113 "Negative"
114 };
115
116 const char * hydroTitles[] =
117 {
118 "Kyte & Doolittle hydropathy parameters",
119 "OHM Hydropathy parameters (Sweet & Eisenberg)",
120 "Consensus parameters (Eisenberg et al)"
121 };
122
123 embInit("pepinfo", argc, argv);
124 ajGraphicsSetPagesize(960, 960);
125
126 ajHistogramSetMark(NOY);
127
128 /* Get parameters */
129 inseq = ajAcdGetSeq("sequence");
130 seq_begin = ajSeqGetBegin(inseq);
131 seq_end = ajSeqGetEnd(inseq);
132 seq_start = seq_begin - 1;
133 outfile = ajAcdGetOutfile("outfile");
134 hwindow = ajAcdGetInt("hwindow");
135 do_general = ajAcdGetBoolean("generalplot");
136 do_hydropathy = ajAcdGetBoolean("hydropathyplot");
137
138 aa_properties = ajAcdGetDatafile("aaproperties");
139 aa_hydropathy = ajAcdGetDatafile("aahydropathy");
140
141 graphs = ajAcdGetGraphxy("graph");
142
143 do_plot = do_general || do_hydropathy;
144
145 /* Set begin and end position in sequence structure */
146 ajSeqSetRange(inseq, seq_begin, seq_end);
147
148 /* Find out which tables are required in the output */
149 do_seq = ajFalse;
150
151 if(do_hydropathy)
152 numGraphs +=3;
153
154 key = ajStrNew();
155 value = ajStrNew();
156
157 /* if sequence plot required */
158 if(do_seq)
159 ajDebug("sequence plot\n");
160
161
162 /* if general properties plot required */
163 if(do_general)
164 {
165 ajDebug("general plot\n");
166
167 /*initialize properties list*/
168 aa_props = ajListNew();
169 embDataListRead(aa_props, aa_properties);
170
171 /* Get first table from properties list of tables */
172 listIter = ajListIterNewread(aa_props);
173
174 /* calculate plot */
175 for(i = 0; i < 9; i++)
176 {
177 if(!ajListIterDone(listIter))
178 {
179 /* ajalloc new ajint array for storing results */
180 AJCNEW(ival,(seq_end-seq_start));
181 iv[i] = ival;
182 table = ajListIterGet(listIter);
183 for(j = seq_start; j < seq_end; j++)
184 {
185 ajStrAssignSubS(&key, ajSeqGetSeqS(inseq), j, j);
186 value = ajTableFetchS(table, key);
187 if(value != NULL)
188 {
189 if(ajStrToInt(value, ival))
190 ival++;
191 else
192 {
193 ajErr("value is not integer ..%S..\n",
194 value);
195 embExit();
196 }
197 }
198 else
199 {
200 ajErr("At position %d in seq, couldn't find key "
201 "%S", j, key);
202 embExit();
203 }
204 }
205 }
206 else
207 {
208 ajErr("No more tables in list\n");
209 embExit();
210 }
211 }
212
213 ajListIterDel(&listIter);
214
215 /* print out results */
216
217 for(i=0; i<9; i++)
218 {
219 ajFmtPrintS(&tmpa, "%s residues in %s from position %d to %d",
220 propertyTitles[i], ajSeqGetNameC(inseq),
221 seq_begin, seq_end);
222 pepinfo_printIntResults(outfile, inseq, iv[i], ajStrGetPtr(tmpa));
223 }
224
225 /* plot out results */
226
227 hist = ajHistNewG(9, (seq_end - seq_begin+1), graphs);
228 hist->bins = seq_end - seq_begin +1;
229
230 hist->xmin = (float) seq_begin;
231 hist->xmax = (float) seq_end;
232
233 hist->displaytype=HIST_SEPARATE;
234
235 ajFmtPrintS(&tmpa, "Properties of residues in %s from position "
236 "%d to %d", ajSeqGetNameC(inseq),seq_begin, seq_end);
237 ajHistSetTitleC(hist, ajStrGetPtr(tmpa));
238
239 ajHistSetXlabelC(hist, "Residue Number");
240 ajHistSetYlabelC(hist, "");
241
242 for(i=0; i<9; i++)
243 {
244 ajFmtPrintS(&tmpa, "%s residues in %s from position %d to %d",
245 propertyTitles[i], ajSeqGetNameC(inseq), seq_begin,
246 seq_end);
247 ajFmtPrintS(&tmpb, "%s residues", propertyTitles[i]);
248 pepinfo_plotHistInt2(hist, iv[i], i,
249 ajStrGetPtr(tmpa), ajStrGetPtr(tmpb), "");
250 }
251
252
253 ajHistDisplay(hist);
254
255
256 for(i = 0; i < 9; i++)
257 AJFREE(iv[i]);
258
259 /*delete hist object*/
260 ajHistDel(&hist);
261 }
262
263 /* if hydropathy plot required */
264 if(do_hydropathy)
265 {
266 ajDebug("hydropathy plot\n");
267
268 if(numGraphs)
269 ajGraphxySetflagOverlay(graphs, ajFalse);
270
271 /* get data from amino acid properties */
272 aa_hydro = ajListNew();
273 embDataListRead(aa_hydro, aa_hydropathy);
274
275 /* Get first table from properties list */
276 listIter = ajListIterNewread(aa_hydro);
277
278 /* calculate plot */
279 for(i=0; i < 3; i++)
280 {
281 /* make sure we have another table from the list to calculate */
282 if(ajListIterDone(listIter))
283 {
284 ajErr("No more tables in list\n");
285 embExit();
286 }
287
288 /* Get next table of parameters */
289 table = ajListIterGet(listIter);
290
291 win_mid = (hwindow / 2);
292
293 /* get array to store result */
294 AJCNEW(pfloat, (seq_end - seq_start));
295 pf[i] = pfloat;
296
297 /* Fill in 0.00 for seq begin to win_mid */
298 for(j=0,cnt=0;j<win_mid; j++)
299 pfloat[cnt++]=0.0;
300
301 /* start loop */
302 for(j = seq_start; j<=(seq_end-hwindow); j++)
303 {
304 total = 0.00;
305 for(k=0; k < hwindow; k++)
306 {
307 ajStrAssignSubS(&key, ajSeqGetSeqS(inseq), (j+k), (j+k));
308 value = ajTableFetchS(table, key);
309 if(value == NULL)
310 {
311 ajErr("At position %d in seq, couldn't find key %s",
312 k, ajStrGetPtr(key));
313 embExit();
314 }
315
316 if(!ajStrIsFloat(value))
317 {
318 ajErr("value is not float ..%s..",
319 ajStrGetPtr(value));
320 embExit();
321 }
322 ajStrToFloat(value, &num);
323 total +=num;
324 }
325 pfloat[cnt++] = total / hwindow;
326 }
327
328 /* fill in value of 0 for end of sequence */
329 for(j = win_mid+1; j<hwindow; j++)
330 pfloat[cnt++] = 0.00;
331 }
332
333 ajListIterDel(&listIter);
334
335 /* Print out results */
336
337 for(i=0; i<3; i++)
338 {
339 ajFmtPrintS(&tmpa, "Results from %s", hydroTitles[i]);
340 pepinfo_printFloatResults(outfile, inseq, pf[i],
341 ajStrGetPtr(tmpa));
342 }
343
344 /*Plot results*/
345 for(i=0; i<3; i++)
346 {
347 ajFmtPrintS(&tmpa,
348 "Hydropathy plot of residues %d to %d of sequence "
349 "%s using %s",seq_begin, seq_end, ajSeqGetNameC(inseq),
350 hydroTitles[i]);
351 pepinfo_plotGraph2Float(graphs, pf[i], ajStrGetPtr(tmpa),
352 "Residue Number", "Hydropathy value");
353 }
354
355 for(i=0; i<3; i++)
356 AJFREE(pf[i]);
357 }
358
359 if(numGraphs)
360 {
361 if(do_general || do_seq)
362 ajGraphNewpage(graphs, ajFalse);
363
364 ajGraphicsSetCharscale(0.50);
365 ajGraphSetTitleC(graphs,"Pepinfo");
366
367 ajGraphxyDisplay(graphs,AJTRUE);
368 }
369
370 if(do_plot)
371 ajGraphicsClose();
372
373 ajSeqDel(&inseq);
374 ajFileClose(&outfile);
375 ajFileClose(&aa_properties);
376 ajFileClose(&aa_hydropathy);
377 ajGraphxyDel(&graphs);
378 ajHistDel(&hist);
379
380 /* Delete Data tables*/
381 embDataListDel(&aa_props);
382 embDataListDel(&aa_hydro);
383
384 ajStrDel(&tmpa);
385 ajStrDel(&tmpb);
386 ajStrDel(&key);
387
388 embExit();
389
390 return 0;
391 }
392
393
394
395
396 /* @funcstatic pepinfo_printIntResults ****************************************
397 **
398 ** prints out a resultsList. Very basic at the moment, really just used to
399 ** prove I have the results and they are correct. There are several of these
400 ** for each primitive data type (as I write them). So far have
401 ** printIntResults and printFloatResults. They are public routines
402 **
403 ** @param [u] outfile [AjPFile] file to output to.
404 ** @param [r] seq [const AjPSeq] Sequence
405 ** @param [r] results [const ajint*] array of reuslts.
406 ** @param [r] header [const char*] header line
407 ** @return [void]
408 ** @@
409 ******************************************************************************/
410
pepinfo_printIntResults(AjPFile outfile,const AjPSeq seq,const ajint * results,const char * header)411 static void pepinfo_printIntResults(AjPFile outfile, const AjPSeq seq,
412 const ajint* results, const char * header)
413 {
414 ajint i;
415 AjPStr aa;
416
417 aa = ajStrNew();
418 ajFmtPrintF(outfile, "Printing out %s\n\n", header);
419 ajFmtPrintF(outfile, "Position Residue\t\t\tResult\n");
420 for(i = seq_start; i<seq_end; i++)
421 {
422 ajStrAssignSubS(&aa, ajSeqGetSeqS(seq), i, i);
423 ajFmtPrintF(outfile, " %5d%8s%32d\n", (i+1),
424 ajStrGetPtr(aa), *results++);
425 }
426
427 ajFmtPrintF(outfile, "\n\n\n");
428
429 ajStrDel(&aa);
430
431 return;
432 }
433
434
435
436
437 /* @funcstatic pepinfo_printFloatResults **************************************
438 **
439 ** Routine to print out Float results data
440 **
441 ** @param [u] outfile [AjPFile] file to output to.
442 ** @param [r] seq [const AjPSeq] Sequence
443 ** @param [r] results [const float*] array of results.
444 ** @param [r] header [const char*] header line
445 ** @return [void]
446 ** @@
447 ******************************************************************************/
448
pepinfo_printFloatResults(AjPFile outfile,const AjPSeq seq,const float * results,const char * header)449 static void pepinfo_printFloatResults(AjPFile outfile, const AjPSeq seq,
450 const float * results,
451 const char * header)
452 {
453 ajint i;
454 AjPStr aa;
455
456 aa = ajStrNew();
457 ajFmtPrintF(outfile, "Printing out %s\n\n", header);
458 ajFmtPrintF(outfile, "Position Residue\t\t\tResult\n");
459 for(i = seq_start; i<seq_end; i++)
460 {
461 ajStrAssignSubS(&aa, ajSeqGetSeqS(seq), i, i);
462 ajFmtPrintF(outfile, "%5d%8s%32.3f\n", (i+1),
463 ajStrGetPtr(aa), *results++);
464 }
465 ajFmtPrintF(outfile, "\n\n\n");
466
467 ajStrDel(&aa);
468
469 return;
470 }
471
472
473
474
475 /* @funcstatic pepinfo_plotGraph2Float ****************************************
476 **
477 ** Create and add graph from set of results to graph set.
478 **
479 ** @param [u] graphs [AjPGraph] Graphs set to add new graph to.
480 ** @param [r] results [const float*] array of results.
481 ** @param [r] title_text [const char*] title for graph
482 ** @param [r] xtext [const char*] x label.
483 ** @param [r] ytext [const char*] y label.
484 ** @return [void]
485 **
486 ******************************************************************************/
487
pepinfo_plotGraph2Float(AjPGraph graphs,const float * results,const char * title_text,const char * xtext,const char * ytext)488 static void pepinfo_plotGraph2Float(AjPGraph graphs,
489 const float * results,
490 const char * title_text,
491 const char * xtext, const char * ytext)
492 {
493
494 AjPGraphdata plot;
495
496 ajint npts = 0;
497
498 float ymin = 64000.;
499 float ymax = -64000.;
500
501 npts = seq_end - seq_start;
502
503 ajGraphicsCalcRange(results,npts,&ymin,&ymax);
504
505 /*
506 ** initialise plot, the number of points will be the length of the data
507 ** in the results structure
508 */
509 plot = ajGraphdataNewI(npts);
510
511 /*Set up rest of plot information*/
512 ajGraphdataSetTitleC(plot, title_text);
513 ajGraphdataSetXlabelC(plot, xtext);
514 ajGraphdataSetYlabelC(plot, ytext);
515 ajGraphdataSetMinmax(plot,(float)1,(float)npts,ymin,ymax);
516 ajGraphdataSetTruescale(plot,(float)1,(float)npts,ymin,ymax);
517 ajGraphdataSetTypeC(plot,"2D Plot");
518
519 ajGraphdataCalcXY(plot, npts, (float)seq_begin, 1.0, results);
520 ajGraphDataAdd(graphs, plot);
521
522 return;
523 }
524
525
526
527
528 /* @funcstatic pepinfo_plotHistInt2 *******************************************
529 **
530 ** Add a histogram data to the set set.
531 **
532 ** @param [u] hist [AjPHist] Histogram set to add new set to.
533 ** @param [r] results [const ajint*] array of results.
534 ** @param [r] hist_num [ajint] the number of the histogram set.
535 ** @param [r] header [const char*] title.
536 ** @param [r] xtext [const char*] x label.
537 ** @param [r] ytext [const char*] y label.
538 ** @return [void]
539 **
540 ******************************************************************************/
541
pepinfo_plotHistInt2(AjPHist hist,const ajint * results,ajint hist_num,const char * header,const char * xtext,const char * ytext)542 static void pepinfo_plotHistInt2(AjPHist hist,
543 const ajint * results,
544 ajint hist_num, const char * header,
545 const char * xtext, const char * ytext)
546 {
547 ajint npts;
548 ajint i;
549
550 float *farray;
551
552 npts = seq_end - seq_start;
553
554 AJCNEW(farray, npts);
555 for(i=0; i<npts; i++)
556 farray[i] = (float) results[i];
557
558 ajHistDataCopy(hist, hist_num, farray);
559
560 ajHistSetmultiTitleC(hist, hist_num, header);
561 ajHistSetmultiXlabelC(hist, hist_num, xtext);
562 ajHistSetmultiYlabelC(hist, hist_num, ytext);
563
564 AJFREE(farray);
565
566 return;
567 }
568