1 /* @source prettyplot application
2 **
3 ** Displays and plots sequence alignments and consensus for PLOTTERS.
4 ** @author Copyright (C) Ian Longden (il@sanger.ac.uk)
5 ** @@
6 **
7 ** Replaces program "prettyplot" (EGCG)
8 **
9 ** options.
10 **
11 ** -ccolours Colour residues by their consensus value. (TRUE)
12 ** -cidentity Colour to display identical matches. (RED)
13 ** -csimilarity Colour to display similar matches. (GREEN)
14 ** -cother Colour to display other matches. (BLACK)
15 **
16 ** -docolour Colour residues by table oily, amide, basic etc. (FALSE)
17 **
18 ** -shade Colour residues by shades. (BLPW)
19 ** B-Black L-Brown P-Wheat W-White
20 **
21 ** -pair values to represent identical similar related (1.5,1.0,0.5)
22 **
23 ** -identity Only match those which are identical in all sequences.
24 **
25 ** -box Display prettybox.
26 **
27 **
28 ** -consensus Display the consensus.
29 **
30 ** -name Display the sequence names.
31 **
32 ** -number Display the residue number at the end of each sequence.
33 **
34 ** -maxnamelen Margin size for the sequence name.
35 **
36 ** -plurality plurality check value used in consensus. (totweight/2)
37 **
38 ** -collision Allow collisions.
39 **
40 ** -portrait Display as a portrait (default landscape).
41 **
42 ** -datafile The data file holding the matrix comparison table.
43 **
44 ** -showscore Obsolete debug variable:
45 ** Print out the scores for a residue number.
46 **
47 ** -alternative 3 other checks for collisions.
48 **
49 **
50 ** This program is free software; you can redistribute it and/or
51 ** modify it under the terms of the GNU General Public License
52 ** as published by the Free Software Foundation; either version 2
53 ** of the License, or (at your option) any later version.
54 **
55 ** This program is distributed in the hope that it will be useful,
56 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
57 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
58 ** GNU General Public License for more details.
59 **
60 ** You should have received a copy of the GNU General Public License
61 ** along with this program; if not, write to the Free Software
62 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
63 ******************************************************************************/
64
65 #include "emboss.h"
66 #include <limits.h>
67
68 #define BOXTOP 0x0001
69 #define BOXBOT 0x0002
70 #define BOXLEF 0x0004
71 #define BOXRIG 0x0008
72 #define BOXCOLOURED 0x0010
73
74 const char **seqcharptr;
75 ajint **seqcolptr;
76 ajint **seqboxptr = NULL;
77 ajint *seqcount = NULL;
78 ajint charlen;
79 AjBool shownames;
80 AjBool shownumbers;
81 AjPSeqset seqset;
82 AjPStr *seqnames;
83 ajint numgaps;
84 char *constr = NULL;
85
86
87
88
89 static ajint prettyplot_calcseqperpage(float yincr,float y,AjBool consensus);
90 static void prettyplot_fillinboxes(ajint length, ajint numseq,
91 ajint start, ajint end,
92 ajint seqstart,ajint seqend,
93 ajint numres, ajint resbreak,
94 AjBool boxit, AjBool boxcol,
95 AjBool consensus,
96 float ystart, float yincr,
97 const AjPSeqCvt cvt);
98
99
100
101
102 /* @prog prettyplot ***********************************************************
103 **
104 ** Displays aligned sequences, with colouring and boxing
105 **
106 ******************************************************************************/
107
main(int argc,char ** argv)108 int main(int argc, char **argv)
109 {
110 ajint i;
111 ajint numseq;
112 ajint j = 0;
113 ajint numres;
114 ajint count;
115 ajint k;
116 ajint kmax;
117 float defheight = 0.0;
118 float currentscale = 0.0;
119 AjPStr shade = NULL;
120 AjPFloat pair = NULL;
121 AjPGraph graph = NULL;
122 AjPMatrix cmpmatrix = NULL;
123 AjPSeqCvt cvt = NULL;
124 AjPStr matcodes = NULL;
125 AjBool consensus;
126 AjBool colourbyconsensus;
127 AjBool colourbyresidues;
128 AjBool colourbyshade = AJFALSE;
129 AjBool boxit;
130 AjBool boxcol;
131 AjBool portrait;
132 AjBool collision;
133 ajint identity;
134 AjBool listoptions;
135 ajint alternative;
136 AjPStr altstr = NULL;
137 AjPStr sidentity = NULL;
138 AjPStr ssimilarity = NULL;
139 AjPStr sother = NULL;
140 AjPStr sboxcolval = NULL;
141 AjPStr options = NULL;
142 /* ajint showscore = 0; */
143 ajint iboxcolval = 0;
144 ajint cidentity = RED;
145 ajint csimilarity = GREEN;
146 ajint cother = BLACK;
147 float fxp;
148 float fyp;
149 float yincr;
150 float y;
151 ajint ixlen;
152 ajint iylen;
153 ajint ixoff;
154 ajint iyoff;
155 char res[2] = " ";
156
157 float *score = 0;
158 float scoremax = 0;
159
160 float *identical = NULL;
161 ajint identicalmaxindex;
162 float *matching = NULL;
163 ajint matchingmaxindex;
164
165 float *colcheck = NULL;
166
167 ajint **matrix;
168 ajint m1 = 0;
169 ajint m2 = 0;
170 ajint ms = 0;
171 ajint highindex = 0;
172 ajint myindex;
173 ajint *previous = 0;
174 AjBool iscons = ajFalse;
175 ajint currentstate = 0;
176 ajint oldfg = 0;
177 float fold = 0.0;
178 ajint *colmat = 0;
179 ajint *shadecolour = 0;
180 /* float identthresh = 1.5; */
181 /* float simthresh = 1.0; */
182 /* float relthresh = 0.5; */
183 float part = 0.0;
184 const char *cptr;
185 ajint resbreak;
186 ajint blocksperline;
187 float ratio;
188 float fplural;
189 float ystart;
190 float xmin;
191 float xmax;
192 float xmid;
193 AjPTime ajtime;
194 ajint gapcount = 0;
195 ajint countforgap = 0;
196 ajint boxindex;
197 float max;
198 ajint matsize;
199 ajint seqperpage = 0;
200 ajint startseq;
201 ajint endseq;
202 ajint newILend = 0;
203 ajint newILstart;
204 void *freeptr;
205 ajint itmp;
206
207 embInit("prettyplot", argc, argv);
208
209 seqset = ajAcdGetSeqset("sequences");
210 numres = ajAcdGetInt("residuesperline");
211 blocksperline = ajAcdGetInt("blocksperline");
212 resbreak = (numres+blocksperline-1)/blocksperline;
213
214 ajSeqsetFill(seqset); /* Pads sequence set with gap characters */
215 numseq = ajSeqsetGetSize(seqset);
216
217 graph = ajAcdGetGraph("graph");
218 colourbyconsensus = ajAcdGetBoolean("ccolours");
219 colourbyresidues = ajAcdGetBoolean("docolour");
220 shade = ajAcdGetString("shade");
221 pair = ajAcdGetArray("pair");
222 identity = ajAcdGetInt("identity");
223 boxit = ajAcdGetBoolean("doboxes");
224
225 ajtime = ajTimeNewTodayFmt("daytime");
226
227 ajSeqsetTrim(seqset);
228 /* offset = ajSeqsetGetOffset(seqset); Unused */
229
230 ajGraphAppendTitleS(graph, ajSeqsetGetUsa(seqset));
231
232 if(boxit)
233 {
234 AJCNEW(seqboxptr, numseq);
235 for(i=0;i<numseq;i++)
236 AJCNEW(seqboxptr[i], ajSeqsetGetLen(seqset));
237 }
238 boxcol = ajAcdGetBoolean("boxcol");
239 sboxcolval = ajAcdGetString("boxuse");
240
241 if(boxcol)
242 {
243 iboxcolval = ajGraphicsCheckColourS(sboxcolval);
244 if(iboxcolval == -1)
245 iboxcolval = GREY;
246 }
247
248 consensus = ajAcdGetBoolean("consensus");
249 if(consensus)
250 {
251 AJCNEW(constr, ajSeqsetGetLen(seqset)+1);
252 constr[0] = '\0';
253 }
254 shownames = ajAcdGetBoolean("name");
255 shownumbers = ajAcdGetBoolean("number");
256 charlen = ajAcdGetInt("maxnamelen");
257 ratio = ajAcdGetFloat("ratio");
258 portrait = ajAcdGetBoolean("portrait");
259 collision = ajAcdGetBoolean("collision");
260 listoptions = ajAcdGetBoolean("listoptions");
261 altstr = ajAcdGetListSingle("alternative");
262 cmpmatrix = ajAcdGetMatrix("matrixfile");
263
264 fplural = ratio * ajSeqsetGetTotweight(seqset);
265
266 ajStrToInt(altstr, &alternative);
267
268 matrix = ajMatrixGetMatrix(cmpmatrix);
269 cvt = ajMatrixGetCvt(cmpmatrix);
270 matsize = ajMatrixGetSize(cmpmatrix);
271
272 AJCNEW(identical,matsize);
273 AJCNEW(matching,matsize);
274 AJCNEW(colcheck,matsize);
275
276 numgaps = numres/resbreak;
277 numgaps--;
278
279 if(portrait)
280 {
281 ajGraphicsSetPortrait(1);
282 ystart = (float) 75.0;
283 }
284 else
285 ystart = (float) 75.0;
286
287 /* pair is an array of three non-negative floats */
288
289 /* identthresh = ajFloatGet(pair,0); Unused */
290 /* simthresh = ajFloatGet(pair,1); Unused */
291 /* relthresh = ajFloatGet(pair,2); Unused */
292
293 /*
294 ** shade is a formatted 4-character string. Characters BLPW only.
295 ** controlled by a pattern in ACD.
296 */
297
298 if(ajStrGetLen(shade))
299 {
300 AJCNEW(shadecolour,4);
301 cptr = ajStrGetPtr(shade);
302 for(i=0;i<4;i++){
303 if(cptr[i]== 'B' || cptr[i]== 'b')
304 shadecolour[i] = BLACK;
305 else if(cptr[i]== 'L' || cptr[i]== 'l')
306 shadecolour[i] = BROWN;
307 else if(cptr[i]== 'P' || cptr[i]== 'p')
308 shadecolour[i] = WHEAT;
309 else if(cptr[i]== 'W' || cptr[i]== 'w')
310 shadecolour[i] = WHITE;
311 }
312
313 colourbyconsensus = colourbyresidues = ajFalse;
314 colourbyshade = ajTrue;
315 }
316
317 /*
318 ** we can colour by consensus or residue but not both
319 ** if we have to choose, use the consensus
320 */
321
322 if(colourbyconsensus && colourbyresidues)
323 colourbyconsensus = AJFALSE;
324
325 sidentity = ajAcdGetString("cidentity");
326 ssimilarity = ajAcdGetString("csimilarity");
327 sother = ajAcdGetString("cother");
328
329 if(colourbyconsensus)
330 {
331 cidentity = ajGraphicsCheckColourS(sidentity);
332 if(cidentity == -1)
333 cidentity = RED;
334
335 csimilarity = ajGraphicsCheckColourS(ssimilarity);
336 if(csimilarity == -1)
337 csimilarity = GREEN;
338
339
340 cother = ajGraphicsCheckColourS(sother);
341 if(cother == -1)
342 cother = BLACK;
343
344 }
345 else if(colourbyresidues)
346 {
347 matcodes = ajMatrixGetCodes(cmpmatrix);
348 if(ajSeqsetIsProt(seqset))
349 colmat = ajGraphicsBasecolourNewProt(matcodes);
350 else
351 colmat = ajGraphicsBasecolourNewNuc(matcodes);
352 }
353
354
355 /* output the options used as the subtitle for the bottom of the graph */
356 if(listoptions)
357 {
358 ajStrAssignC(&options,"");
359 ajFmtPrintAppS(&options,"-plurality %.1f",fplural);
360
361 if(collision)
362 ajStrAppendC(&options," -collision");
363 else
364 ajStrAppendC(&options," -nocollision");
365
366 if(boxit)
367 ajStrAppendC(&options," -box");
368 else
369 ajStrAppendC(&options," -nobox");
370
371 if(boxcol)
372 ajStrAppendC(&options," -boxcol");
373 else
374 ajStrAppendC(&options," -noboxcol");
375
376 if(colourbyconsensus)
377 ajStrAppendC(&options," -colbyconsensus");
378 else if(colourbyresidues)
379 ajStrAppendC(&options," -colbyresidues");
380 else if(colourbyshade)
381 ajStrAppendC(&options," -colbyshade");
382 else
383 ajStrAppendC(&options," -nocolour");
384
385 if(alternative==2)
386 ajStrAppendC(&options," -alt 2");
387 else if(alternative==1)
388 ajStrAppendC(&options," -alt 1");
389 else if(alternative==3)
390 ajStrAppendC(&options," -alt 3");
391 }
392
393
394 AJCNEW(seqcolptr, numseq);
395 for(i=0;i<numseq;i++)
396 AJCNEW(seqcolptr[i], ajSeqsetGetLen(seqset));
397
398 AJCNEW(seqcharptr, numseq);
399 AJCNEW(seqnames, numseq);
400 AJCNEW(score, numseq);
401 AJCNEW(previous, numseq);
402 AJCNEW(seqcount, numseq);
403
404 for(i=0;i<numseq;i++)
405 {
406 ajSeqsetFmtUpper(seqset);
407 seqcharptr[i] = ajSeqsetGetseqSeqC(seqset, i);
408 seqnames[i] = 0;
409 ajStrAppendS(&seqnames[i],ajSeqsetGetseqNameS(seqset, i));
410 ajStrTruncateLen(&seqnames[i],charlen);
411 previous[i] = 0;
412 seqcount[i] = 0;
413 }
414
415 /*
416 ** user will pass the number of residues to fit a page
417 ** therefore we now need to calculate the size of the chars
418 ** based on this and get the new char width.
419 ** 'charlen' maximum characters for the name (truncated above)
420 */
421
422 ajGraphicsGetCharsize(&defheight,¤tscale);
423
424 xmin = -charlen - (float)2.0;
425 xmax = (float)numres+(float)11.0+(float)(numres/resbreak);
426 xmid = (xmax + xmin)/(float)2.0;
427
428 ajGraphOpenWin(graph, xmin, xmax,
429 (float)0.0, ystart+(float)1.0);
430
431 ajGraphGetParamsPage(graph, &fxp,&fyp,&ixlen,&iylen,&ixoff,&iyoff);
432
433 if(portrait)
434 {
435 itmp = ixlen;
436 ixlen = iylen;
437 iylen = itmp;
438 }
439
440 ajGraphicsGetCharsize(&defheight,¤tscale);
441
442 ajGraphicsSetCharscale(((float)ixlen/((float)(numres+charlen+1)*
443 (currentscale * (float) 1.5)))/
444 currentscale);
445
446 /* ajGraphicsSetCharscale(((float)ixlen/((float)(numres+charlen)*
447 (currentscale+(float)1.0)))/
448 currentscale); */
449
450 ajGraphicsGetCharsize(&defheight,¤tscale);
451
452 yincr = (currentscale + (float)3.0)*(float)0.3;
453
454 /*
455 ** If we have titles (now the standard graph title and subtitle and footer)
456 ** leave 7 rows of space for them
457 */
458 y=ystart-(float)7.0;
459
460 if(ajStrGetLen(options))
461 {
462 fold = ajGraphicsSetCharscale(1.0);
463 ajGraphicsDrawposTextAtmid(xmid,2.0,
464 ajStrGetPtr(options));
465 ajGraphicsSetCharscale(fold);
466 }
467
468 /* if sequences per page not set then calculate it */
469
470 if(!seqperpage)
471 {
472 seqperpage = prettyplot_calcseqperpage(yincr,y,consensus);
473 if(seqperpage>numseq)
474 seqperpage=numseq;
475 }
476
477 count = 0;
478
479 /*
480 ** for boxes we need to set a foreground colour for the box lines
481 ** and save the current foreground colour
482 */
483 if(boxit && boxcol)
484 oldfg = ajGraphicsSetFgcolour(iboxcolval);
485
486 /*
487 ** step through each residue position
488 */
489
490 kmax = ajSeqsetGetLen(seqset) - 1;
491 for(k=0; k<= kmax; k++)
492 {
493 /* reset column score array */
494 for(i=0;i<numseq;i++)
495 score[i] = 0.0;
496
497 /* reset matrix character testing arrays */
498 for(i=0;i<matsize;i++)
499 {
500 identical[i] = 0.0;
501 matching[i] = 0.0;
502 colcheck[i] = 0.0;
503 }
504
505 /* generate a score for this residue in each sequence */
506 for(i=0;i<numseq;i++)
507 {
508 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
509 for(j=0;j<numseq;j++)
510 {
511 m2 = ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]);
512 if(m1 && m2)
513 score[i] += (float)matrix[m1][m2]*
514 ajSeqsetGetseqWeight(seqset, j);
515 }
516 if(m1)
517 identical[m1] += ajSeqsetGetseqWeight(seqset, i);
518 }
519
520 /* find the highest score */
521 highindex = -1;
522 scoremax = INT_MIN;
523 /*ajDebug("Scores at position %d:\n", k);*/
524
525 for(i=0;i<numseq;i++)
526 {
527 /*ajDebug(" seq %d: '%c' %f\n",i,seqcharptr[i][k],score[i]);*/
528
529 if(score[i] > scoremax)
530 {
531 scoremax = score[i];
532 highindex = i;
533 }
534 }
535 for(i=0;i<numseq;i++)
536 {
537 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
538
539 if(!matching[m1])
540 {
541 for(j=0;j<numseq;j++)
542 {
543 m2 = ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]);
544 if(m1 && m2 && matrix[m1][m2] > 0)
545 matching[m1] += ajSeqsetGetseqWeight(seqset, j);
546 }
547 }
548 }
549
550 /* find highs for matching and identical */
551 matchingmaxindex = 0;
552 identicalmaxindex = 0;
553 for(i=0;i<numseq;i++)
554 {
555 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
556 if(identical[m1] > identical[identicalmaxindex])
557 identicalmaxindex = m1;
558 }
559 for(i=0;i<numseq;i++)
560 {
561 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
562 if(matching[m1] > matching[matchingmaxindex])
563 matchingmaxindex = m1;
564 else if(matching[m1] == matching[matchingmaxindex])
565 {
566 if(identical[m1] > identical[matchingmaxindex])
567 matchingmaxindex= m1;
568 }
569 }
570
571 iscons = ajFalse;
572 boxindex = -1;
573 max = -3;
574
575 ajDebug("k:%2d highindex:%2d matching:%4.2f\n",
576 k, highindex,
577 matching[ajSeqcvtGetCodeK(cvt, seqcharptr[highindex][k])]);
578 if(highindex != -1 &&
579 matching[ajSeqcvtGetCodeK(cvt, seqcharptr[highindex][k])] >= fplural)
580 {
581 iscons = ajTrue;
582 boxindex = highindex;
583 }
584 else
585 {
586 for(i=0;i<numseq;i++)
587 {
588 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
589 if(matching[m1] > max)
590 {
591 max = matching[m1];
592 highindex = i;
593 }
594 else if(matching[m1] == max)
595 {
596 if(identical[m1] >
597 identical[ajSeqcvtGetCodeK(cvt,
598 seqcharptr[highindex][k])] )
599 {
600 max = matching[m1];
601 highindex = i;
602 }
603 }
604 }
605
606 if(matching[ajSeqcvtGetCodeK(cvt,
607 seqcharptr[highindex][k])] >= fplural)
608 {
609 iscons = ajTrue;
610 boxindex = highindex;
611 }
612 }
613
614
615 if(iscons)
616 {
617 if(!collision)
618 {
619 /* check for collisions */
620 if(alternative == 1)
621 {
622 /* check to see if this is unique for collisions */
623 for(i=0;i<numseq;i++)
624 {
625 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
626 if(identical[m1] >= identical[identicalmaxindex] &&
627 m1 != identicalmaxindex)
628 iscons = ajFalse;
629 }
630
631 /*ajDebug("after (alt=1) iscons: %B",iscons);*/
632 }
633
634 else if(alternative == 2)
635 {
636 for(i=0;i<numseq;i++)
637 {
638 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
639
640 if((matching[m1] >= matching[matchingmaxindex] &&
641 m1 != matchingmaxindex &&
642 matrix[m1][matchingmaxindex] < 0.1)||
643 (identical[m1] >= identical[matchingmaxindex]
644 && m1 != matchingmaxindex))
645 iscons = ajFalse;
646 }
647 }
648 else if(alternative == 3)
649 {
650 /*
651 ** to do this check one is NOT in consensus to see if
652 ** another score of fplural has been found
653 */
654 ms = ajSeqcvtGetCodeK(cvt, seqcharptr[highindex][k]);
655
656 for(i=0;i<numseq;i++)
657 {
658 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
659 if(ms != m1 && colcheck[m1] == 0.0)
660 /* NOT in the current consensus */
661 for(j=0;j<numseq;j++)
662 {
663 m2 = ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]);
664 if( matrix[ms][m2] < 0.1)
665 {
666 /* NOT in the current consensus */
667 if( matrix[m1][m2] > 0.1)
668 colcheck[m1] +=
669 ajSeqsetGetseqWeight(seqset,
670 j);
671 }
672 }
673 }
674
675 for(i=0;i<numseq;i++)
676 {
677 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
678 /* if any other matches then we have a collision */
679 if(colcheck[m1] >= fplural)
680 iscons = ajFalse;
681 }
682
683 /*ajDebug("after alt=2 iscons: %B", iscons);*/
684 }
685 else
686 {
687 for(i=0;i<numseq;i++)
688 {
689 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
690 if((matching[m1] >= matching[matchingmaxindex] &&
691 m1 != matchingmaxindex &&
692 matrix[m1][matchingmaxindex] < 0.1))
693 iscons = ajFalse;
694 if(identical[m1] >= identical[matchingmaxindex] &&
695 m1 != matchingmaxindex &&
696 matrix[m1][matchingmaxindex] > 0.1)
697 iscons = ajFalse;
698 }
699
700 if(!iscons)
701 { /* matches failed try identicals */
702 if(identical[identicalmaxindex] >= fplural)
703 {
704 iscons = ajTrue;
705 /*
706 ** if nothing has an equal or higher match that
707 ** does not match highest then false
708 */
709 for(i=0;i<numseq;i++)
710 {
711 m1 = ajSeqcvtGetCodeK(cvt, seqcharptr[i][k]);
712 if(identical[m1] >=
713 identical[identicalmaxindex] &&
714 m1 != identicalmaxindex)
715 iscons = ajFalse;
716 else if(matching[m1] >=
717 matching[identicalmaxindex] &&
718 matrix[m1][matchingmaxindex] <= 0.0)
719 iscons = ajFalse;
720 else if(m1 == identicalmaxindex)
721 j = i;
722 }
723
724 if(iscons)
725 highindex = j;
726 }
727 }
728
729 }
730 }
731
732 if(identity)
733 {
734 j = 0;
735 for(i=0;i<numseq;i++)
736 if(seqcharptr[highindex][k] == seqcharptr[i][k])
737 j++;
738
739 if(j<identity)
740 iscons = ajFalse;
741 }
742 }
743
744 /*
745 ** Done a full line of residues
746 ** Boxes have been defined up to this point
747 */
748 if(count >= numres )
749 {
750 /* check y position for next set */
751 y=y-(yincr*((float)numseq+(float)2.0+((float)consensus*(float)2)));
752 if(y<yincr*((float)numseq+(float)2.0+((float)consensus*(float)2)))
753 {
754 /* full page - print it */
755 y=ystart-(float)6.0;
756
757 startseq = 0;
758 endseq = seqperpage;
759 newILstart = newILend;
760 newILend = k;
761 while(startseq < numseq)
762 {
763 /* AJB */
764 /*if(startseq != 0)
765 ajGraphNewpage(graph, AJFALSE);*/
766
767 /*ajDebug("Inner loop: startseq: %d numseq: %d endseq: %d\n",
768 startseq, numseq, endseq);*/
769 if(endseq>numseq)
770 endseq=numseq;
771 prettyplot_fillinboxes(numseq,ajSeqsetGetLen(seqset),
772 startseq,endseq,
773 newILstart,newILend,
774 numres,resbreak,
775 boxit,boxcol,consensus,
776 ystart,yincr,cvt);
777 startseq = endseq;
778 endseq += seqperpage;
779 ajGraphNewpage(graph, AJFALSE);
780 }
781 }
782
783 count = 0;
784 gapcount = 0;
785 }
786
787 count++;
788 countforgap++;
789
790 for(j=0;j<numseq;j++)
791 {
792 /* START OF BOXES */
793
794 if(boxit)
795 {
796 seqboxptr[j][k] = 0;
797 if(boxindex!=-1)
798 {
799 myindex = boxindex;
800 if(matrix[ajSeqcvtGetCodeK(cvt, seqcharptr[j][k])]
801 [ajSeqcvtGetCodeK(cvt, seqcharptr[myindex][k])] > 0)
802 part = 1.0;
803 else
804 {
805 if(identical[ajSeqcvtGetCodeK(cvt,
806 seqcharptr[j][k])] >=
807 fplural)
808 part = 1.0;
809 else
810 part = 0.0;
811 }
812
813 if(previous[j] != part)
814 /* draw vertical line */
815 seqboxptr[j][k] |= BOXLEF;
816
817 if(j==0)
818 {
819 /* special case for horizontal line */
820 if(part)
821 {
822 currentstate = 1;
823 /* draw hori line */
824 seqboxptr[j][k] |= BOXTOP;
825 }
826 else
827 currentstate = 0;
828 }
829 else
830 {
831 /* j != 0 Normal case for horizontal line */
832 if(part != currentstate)
833 {
834 /*draw hori line */
835 seqboxptr[j][k] |= BOXTOP;
836 currentstate = (ajint) part;
837 }
838 }
839
840 if(j== numseq-1 && currentstate)
841 /* draw horiline at bottom */
842 seqboxptr[j][k] |= BOXBOT;
843
844 previous[j] = (ajint) part;
845 }
846 else
847 {
848 part = 0;
849 if(previous[j])
850 {
851 /* draw vertical line */
852 seqboxptr[j][k] |= BOXLEF;
853 }
854 previous[j] = 0;
855 }
856
857 if(count == numres || k == kmax || countforgap >= resbreak )
858 { /* last one on the row or a break*/
859 if(previous[j])
860 {
861 /* draw vertical line */
862 seqboxptr[j][k] |= BOXRIG;
863 }
864 previous[j] = 0;
865 }
866
867 } /* end box */
868
869 if(boxit && boxcol)
870 if(boxindex != -1)
871 {
872 myindex = boxindex;
873 if(matrix[ajSeqcvtGetCodeK(cvt, seqcharptr[j][k])]
874 [ajSeqcvtGetCodeK(cvt, seqcharptr[myindex][k])] > 0
875 || identical[ajSeqcvtGetCodeK(cvt, seqcharptr[j][k])] >=
876 fplural )
877
878 seqboxptr[j][k] |= BOXCOLOURED;
879 }
880
881 /* END OF BOXES */
882
883
884
885
886 if(ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]))
887 res[0] = seqcharptr[j][k];
888 else
889 res[0] = '-';
890
891 if(colourbyconsensus)
892 {
893 part = (float) matrix[ajSeqcvtGetCodeK(cvt, seqcharptr[j][k])]
894 [ajSeqcvtGetCodeK(cvt, seqcharptr[highindex][k])];
895 if(iscons && seqcharptr[highindex][k] == seqcharptr[j][k])
896 seqcolptr[j][k] = cidentity;
897 else if(part > 0.0)
898 seqcolptr[j][k] = csimilarity;
899 else
900 seqcolptr[j][k] = cother;
901 }
902 else if(colourbyresidues)
903 seqcolptr[j][k] = colmat[ajSeqcvtGetCodeK(cvt,
904 seqcharptr[j][k])];
905 else if(iscons && colourbyshade)
906 {
907 part = (float) matrix[ajSeqcvtGetCodeK(cvt, seqcharptr[j][k])]
908 [ajSeqcvtGetCodeK(cvt, seqcharptr[highindex][k])];
909 if(part >= 1.5)
910 seqcolptr[j][k] = shadecolour[0];
911 else if(part >= 1.0)
912 seqcolptr[j][k] = shadecolour[1];
913 else if(part >= 0.5)
914 seqcolptr[j][k] = shadecolour[2];
915 else
916 seqcolptr[j][k] = shadecolour[3];
917 }
918 else if(colourbyshade)
919 seqcolptr[j][k] = shadecolour[3];
920 else
921 seqcolptr[j][k] = BLACK;
922 }
923
924 if(consensus)
925 {
926 if(iscons)
927 res[0] = seqcharptr[highindex][k];
928 else
929 res[0] = '-';
930 strcat(constr,res);
931 }
932
933 if(countforgap >= resbreak)
934 {
935 gapcount++;
936 countforgap=0;
937 }
938 }
939
940
941 startseq = 0;
942 endseq=seqperpage;
943 newILstart = newILend;
944 newILend = k;
945 while(startseq < numseq)
946 {
947 if(startseq)
948 ajGraphNewpage(graph, AJFALSE);
949
950 /*ajDebug("Final loop: startseq: %d numseq: %d endseq: %d\n",
951 startseq, numseq, endseq);*/
952 if(endseq>numseq)
953 endseq = numseq;
954 prettyplot_fillinboxes(numseq,ajSeqsetGetLen(seqset),
955 startseq,endseq,
956 newILstart,newILend,
957 numres,resbreak,
958 boxit,boxcol,consensus,
959 ystart,yincr,cvt);
960 startseq = endseq;
961 endseq += seqperpage;
962 }
963
964
965 ajGraphicsGetCharsize(&defheight,¤tscale);
966
967 if(boxit && boxcol)
968 oldfg = ajGraphicsSetFgcolour(oldfg);
969
970 ajGraphicsCloseWin();
971 ajGraphxyDel(&graph);
972
973 ajStrDel(&sidentity);
974 ajStrDel(&ssimilarity);
975 ajStrDel(&sother);
976 ajStrDel(&options);
977 ajStrDel(&altstr);
978
979 ajStrDel(&matcodes);
980
981 for(i=0;i<numseq;i++)
982 {
983 ajStrDel(&seqnames[i]);
984 AJFREE(seqcolptr[i]);
985 if(seqboxptr)
986 AJFREE(seqboxptr[i]);
987 }
988 AJFREE(seqcolptr);
989 AJFREE(seqboxptr);
990
991 AJFREE(seqnames);
992 AJFREE(score);
993 AJFREE(previous);
994 AJFREE(seqcount);
995
996 AJFREE(colmat);
997 AJFREE(shadecolour);
998
999 freeptr = (void *) seqcharptr;
1000 AJFREE(freeptr);
1001
1002 AJFREE(identical);
1003 AJFREE(matching);
1004 AJFREE(colcheck);
1005
1006 ajSeqsetDel(&seqset);
1007 ajMatrixDel(&cmpmatrix);
1008 ajStrDel(&shade);
1009 ajStrDel(&sboxcolval);
1010 ajStrDel(&sidentity);
1011 ajStrDel(&ssimilarity);
1012 ajStrDel(&sother);
1013 ajFloatDel(&pair);
1014 ajTimeDel(&ajtime);
1015 AJFREE(constr);
1016
1017 embExit();
1018
1019 return 0;
1020 }
1021
1022
1023
1024
1025 /* @funcstatic prettyplot_calcseqperpage **************************************
1026 **
1027 ** Calculate the number of sequences per page.
1028 **
1029 ** @param [r] yincr [float] Undocumented
1030 ** @param [r] y [float] Undocumented
1031 ** @param [r] consensus [AjBool] If true, include a consensus sequence
1032 ** @return [ajint] Number of sequences that fit on one page
1033 ** @@
1034 ******************************************************************************/
1035
1036
prettyplot_calcseqperpage(float yincr,float y,AjBool consensus)1037 static ajint prettyplot_calcseqperpage(float yincr,float y,AjBool consensus)
1038 {
1039 float yep = 1.0;
1040 ajint numallowed = 1;
1041
1042 while(yep>0.0)
1043 {
1044 yep = y-(yincr*((float)numallowed+(float)2.0+((float)consensus*
1045 (float)2)));
1046 numallowed++;
1047 }
1048
1049 return numallowed-1;
1050 }
1051
1052
1053
1054
1055 /* @funcstatic prettyplot_fillinboxes *****************************************
1056 **
1057 ** Plot the page
1058 **
1059 ** @param [r] numseq [ajint] Number of sequences in alignment
1060 ** @param [r] length [ajint] Alignment length
1061 ** @param [r] seqstart [ajint] First sequence to plot
1062 ** @param [r] seqend [ajint] Last sequence to plot
1063 ** @param [r] start [ajint] First residue
1064 ** @param [r] end [ajint] Last residue
1065 ** @param [r] numres [ajint] Number of residues per line
1066 ** @param [r] resbreak [ajint] Number of residues per block
1067 ** @param [r] boxit [AjBool] If true, display box outline
1068 ** @param [r] boxcol [AjBool] If true, colour the background in each box
1069 ** @param [r] consensus [AjBool] If true, include consensus sequence
1070 ** on last line
1071 ** @param [r] ystart [float] Vertical start. 1.0 is the top.
1072 ** @param [r] yincr [float] Vertical increment to next row
1073 ** @param [r] cvt [const AjPSeqCvt] Conversion table from residue code
1074 ** to matrix position
1075 ** @return [void]
1076 ** @@
1077 ******************************************************************************/
1078
prettyplot_fillinboxes(ajint numseq,ajint length,ajint seqstart,ajint seqend,ajint start,ajint end,ajint numres,ajint resbreak,AjBool boxit,AjBool boxcol,AjBool consensus,float ystart,float yincr,const AjPSeqCvt cvt)1079 static void prettyplot_fillinboxes(ajint numseq, ajint length,
1080 ajint seqstart, ajint seqend,
1081 ajint start, ajint end,
1082 ajint numres, ajint resbreak,
1083 AjBool boxit, AjBool boxcol,
1084 AjBool consensus,
1085 float ystart, float yincr,
1086 const AjPSeqCvt cvt)
1087 {
1088 ajint count = 1;
1089 ajint gapcount = 0;
1090 ajint countforgap = 0;
1091 ajint table[16];
1092 ajint i;
1093 ajint j;
1094 ajint k;
1095 ajint w;
1096 ajint oldfg = 0;
1097 ajint oldcol = 0;
1098 ajint l;
1099 float y;
1100 char res[2]=" ";
1101 AjPStr strcon = NULL;
1102 char numberstring[10];
1103 float defcs = 0.;
1104 float curcs = 0.;
1105
1106 /*
1107 ajDebug("fillinboxes numseq: %d length: %d\n",
1108 numseq, length);
1109 ajDebug("fillinboxes start: %d end: %d seqstart: %d seqend; %d\n",
1110 start, end, seqstart, seqend);
1111 ajDebug("fillinboxes numres: %d resbreak: %d\n",
1112 numres, resbreak);
1113 ajDebug("fillinboxes boxit: %B boxcol: %B consensus: %B \n",
1114 boxit, boxcol, consensus);
1115 ajDebug("fillinboxes xmid: %.3f ystart:%.3f yincr: %.3f\n",
1116 xmid, ystart, yincr);
1117 */
1118 ajStrAppendC(&strcon,"Consensus");
1119 ajStrTruncateLen(&strcon,charlen);
1120
1121 if(boxit && boxcol)
1122 {
1123 y = ystart-(float)6.0;
1124 for(k=start; k< end; k++)
1125 {
1126 if(countforgap >= resbreak)
1127 {
1128 gapcount++;
1129 countforgap = 0;
1130 }
1131 if(count >= numres+1 )
1132 {
1133 y = y-(yincr*((float)numseq+(float)2.0+((float)consensus*
1134 (float)2)));
1135 if(y<yincr*((float)numseq+(float)2.0+((float)consensus*
1136 (float)2)))
1137 {
1138 y = ystart-(float)6.0;
1139 }
1140 count = 1;
1141 gapcount = 0;
1142 }
1143 count++;
1144 countforgap++;
1145 /* thiscol = ajGraphicsGetColourFore(); Unused */
1146
1147 for(j=seqstart,l=0;j<seqend;j++,l++)
1148 if(seqboxptr[j][k] & BOXCOLOURED)
1149 {
1150 ajGraphicsDrawposRectFill((float)(count+gapcount-(float)1)+
1151 (float)1.0,
1152 y-(yincr*((float)l+(float)0.5)),
1153 (float)(count+gapcount-(float)1),
1154 y-(yincr*((float)l-(float)0.5)));
1155 }
1156 }
1157 }
1158
1159 oldcol = ajGraphicsSetFgcolour(BLACK);
1160
1161 /* DO THE BACKGROUND OF THE BOXES FIRST */
1162
1163 count = 0;
1164 gapcount = countforgap = 0;
1165 y = ystart-(float)6.0;
1166
1167 ajGraphicsGetCharsize(&defcs,&curcs);
1168
1169 if(shownames)
1170 {
1171 for(i=seqstart,l=0;i<seqend;i++,l++)
1172 ajGraphicsDrawposTextAtstart((float)-charlen,y-(yincr*l),
1173 ajStrGetPtr(seqnames[i]));
1174
1175 if(consensus && (numseq==seqend))
1176 ajGraphicsDrawposTextAtstart(
1177 (float)-charlen,
1178 y-(yincr*((seqend-seqstart)+(float)1)),ajStrGetPtr(strcon));
1179
1180 }
1181
1182 for(k=start; k< end; k++)
1183 {
1184 if(countforgap >= resbreak)
1185 {
1186 gapcount++;
1187 countforgap = 0;
1188 }
1189
1190 if(count >= numres )
1191 {
1192 if(shownumbers)
1193 {
1194 for(j=seqstart,l=0;j<seqend;j++,l++)
1195 {
1196 sprintf(numberstring,"%d",seqcount[j]);
1197 ajGraphicsDrawposTextAtstart(
1198 (float)(count+numgaps)+(float)5.0,
1199 y-(yincr*(float)l),numberstring);
1200 }
1201
1202 if(consensus && (numseq==seqend))
1203 {
1204 sprintf(numberstring,"%d",k);
1205 ajGraphicsDrawposTextAtstart(
1206 (float)(count+numgaps)+(float)5.0,
1207 y-(yincr*((float)l+(float)1)),
1208 numberstring);
1209 }
1210 }
1211
1212 y = y-(yincr*((float)numseq+(float)2.0+((float)consensus*
1213 (float)2)));
1214 if(y<yincr*((float)numseq+(float)2.0+((float)consensus*(float)2)))
1215 {
1216 y = ystart-(float)6.0;
1217 }
1218
1219 count = 0;
1220 gapcount = 0;
1221 if(shownames)
1222 {
1223 for(i=seqstart,l=0;i<seqend;i++,l++)
1224 ajGraphicsDrawposTextAtstart((float)-charlen,y-(yincr*l),
1225 ajStrGetPtr(seqnames[i]));
1226
1227 if(consensus &&(numseq==seqend))
1228 ajGraphicsDrawposTextAtstart(
1229 (float)-charlen,
1230 y-(yincr*((seqend-seqstart)+1)),
1231 ajStrGetPtr(strcon));
1232 }
1233 }
1234 count++;
1235 countforgap++;
1236
1237 if(boxit)
1238 {
1239 for(j=seqstart,l=0; j< seqend; j++,l++)
1240 {
1241 if(ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]))
1242 seqcount[j]++;
1243 if(seqboxptr[j][k] & BOXLEF)
1244 ajGraphicsDrawposLine((float)(count+gapcount),y-
1245 (yincr*((float)l+(float)0.5)),
1246 (float)(count+gapcount),
1247 y-(yincr*((float)l-(float)0.5)));
1248
1249 if(seqboxptr[j][k] & BOXTOP)
1250 ajGraphicsDrawposLine((float)(count+gapcount),y-
1251 (yincr*((float)l-(float)0.5)),
1252 (float)(count+gapcount)+(float)1.0,
1253 y-(yincr*((float)l-(float)0.5)));
1254
1255 if(seqboxptr[j][k] & BOXBOT)
1256 ajGraphicsDrawposLine((float)(count+gapcount),y-
1257 (yincr*((float)l+(float)0.5)),
1258 (float)(count+gapcount)+(float)1.0,
1259 y-(yincr*((float)l+(float)0.5)));
1260
1261 if(seqboxptr[j][k] & BOXRIG)
1262 ajGraphicsDrawposLine((float)(count+gapcount)+(float)1.0,y-
1263 (yincr*((float)l+(float)0.5)),
1264 (float)(count+gapcount)+(float)1.0,
1265 y-(yincr*((float)l-(float)0.5)));
1266 }
1267 }
1268 else if(shownumbers) /* usually set in the boxit block */
1269 {
1270 for(j=seqstart,l=0; j< seqend; j++,l++)
1271 {
1272 if(ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]))
1273 seqcount[j]++;
1274 }
1275 }
1276 if(consensus && (numseq==seqend))
1277 {
1278 res[0] = constr[k];
1279
1280 /* start -> mid */
1281 ajGraphicsDrawposTextAtmid((float) 0.5 + (float)(count+gapcount),
1282 y-(yincr*((seqend-seqstart)+1)),res);
1283 }
1284 }
1285
1286 if(shownumbers)
1287 {
1288 for(j=seqstart,l=0;j<seqend;j++,l++)
1289 {
1290 sprintf(numberstring,"%d",seqcount[j]);
1291 ajGraphicsDrawposTextAtstart((float)(count+numgaps)+
1292 (float)5.0,y-(yincr*l),
1293 numberstring);
1294 }
1295
1296 if(consensus && (numseq==seqend))
1297 {
1298 sprintf(numberstring,"%d",k);
1299 ajGraphicsDrawposTextAtstart((float)(count+numgaps)+(float)5.0,
1300 y-(yincr*(l+(float)1)),
1301 numberstring);
1302 }
1303 }
1304
1305
1306 ajStrDel(&strcon);
1307
1308 for(i=0;i<16;i++)
1309 table[i] = -1;
1310 for(i=0;i<numseq;i++)
1311 for(k=0; k< length; k++)
1312 table[seqcolptr[i][k]] = 1;
1313
1314 /* now display again but once for each colour */
1315
1316 /* for(w=0;w<15;w++)*/
1317 for(w=0;w<16;w++)
1318 {
1319 /* not 16 as we can ignore white on plotters*/
1320 if(table[w] > 0)
1321 {
1322 oldfg = ajGraphicsSetFgcolour(w);
1323 count = 0;
1324 gapcount = countforgap = 0;
1325
1326 y = ystart-(float)6.0;
1327
1328 for(k=start; k< end; k++)
1329 {
1330 if(countforgap >= resbreak)
1331 {
1332 gapcount++;
1333 countforgap=0;
1334 }
1335
1336 if(count >= numres )
1337 {
1338 y=y-(yincr*((float)(seqend-seqstart)+(float)2.0+
1339 ((float)consensus*(float)2)));
1340 count = 0;
1341 gapcount = 0;
1342 }
1343 count++;
1344 countforgap++;
1345
1346 for(j=seqstart,l=0; j< seqend; j++,l++){
1347 if(seqcolptr[j][k] == w)
1348 {
1349 if(ajSeqcvtGetCodeK(cvt, seqcharptr[j][k]))
1350 res[0] = seqcharptr[j][k];
1351 else
1352 res[0] = '-';
1353 /* start -> mid */
1354 ajGraphicsDrawposTextAtmid((float) 0.5 +
1355 (float) (count+gapcount),
1356 y-(yincr*l),res);
1357 }
1358 }
1359 }
1360 }
1361 oldfg = ajGraphicsSetFgcolour(oldfg);
1362 }
1363
1364 oldfg = ajGraphicsSetFgcolour(oldcol);
1365 start = end;
1366
1367 return;
1368 }
1369