1 /* @source dotmatcher application
2 **
3 ** dotmatcher displays a dotplot for two sequences.
4 **
5 ** @author Copyright (C) Ian Longden (il@sanger.ac.uk)
6 ** @modified: Added non-proportional plot. Copyright (C) Alan Bleasby
7 ** @@
8 **
9 ** This program is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU General Public License
11 ** as published by the Free Software Foundation; either version 2
12 ** of the License, or (at your option) any later version.
13 **
14 ** This program is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ** GNU General Public License for more details.
18 **
19 ** You should have received a copy of the GNU General Public License
20 ** along with this program; if not, write to the Free Software
21 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 ******************************************************************************/
23
24 #include "emboss.h"
25
26
27
28
29 static void dotmatcher_pushpoint(AjPList *l, float x1, float y1, float x2,
30 float y2, AjBool stretch);
31
32
33
34
35 /* @datastatic PPoint *********************************************************
36 **
37 ** Dotmatcher point data
38 **
39 ** @alias SPoint
40 ** @alias OPoint
41 **
42 ** @attr x1 [float] x1 coordinate
43 ** @attr y1 [float] y1 coordinate
44 ** @attr x2 [float] x2 coordinate
45 ** @attr y2 [float] y2 coordinate
46 ******************************************************************************/
47
48 typedef struct SPoint
49 {
50 float x1;
51 float y1;
52 float x2;
53 float y2;
54 } OPoint;
55 #define PPoint OPoint*
56
57
58
59
60 /* @prog dotmatcher ***********************************************************
61 **
62 ** Displays a thresholded dotplot of two sequences
63 **
64 ******************************************************************************/
65
main(int argc,char ** argv)66 int main(int argc, char **argv)
67 {
68 AjPList list = NULL;
69 AjPSeq seq;
70 AjPSeq seq2;
71 AjPStr aa0str = 0;
72 AjPStr aa1str = 0;
73 const char *s1;
74 const char *s2;
75 char *strret = NULL;
76 ajuint i;
77 ajuint j;
78 ajuint k;
79 ajint l;
80 ajint abovethresh;
81 ajint total;
82 ajint starti = 0;
83 ajint startj = 0;
84 ajint windowsize;
85 float thresh;
86 AjPGraph graph = NULL;
87 AjPGraph xygraph = NULL;
88 float flen1;
89 float flen2;
90 ajuint len1;
91 ajuint len2;
92
93 AjPTime ajtime = NULL;
94 time_t tim;
95 AjBool boxit=AJTRUE;
96 /* Different ticks as they need to be different for x and y due to
97 length of string being important on x */
98 ajuint acceptableticksx[]=
99 {
100 1,10,50,100,500,1000,1500,10000,
101 500000,1000000,5000000
102 };
103 ajuint acceptableticks[]=
104 {
105 1,10,50,100,200,500,1000,2000,5000,10000,15000,
106 500000,1000000,5000000
107 };
108 ajint numbofticks = 10;
109 float xmargin;
110 float ymargin;
111 float ticklen;
112 float tickgap;
113 float onefifth;
114 float k2;
115 float max;
116 char ptr[10];
117 AjPMatrix matrix = NULL;
118 ajint** sub;
119 AjPSeqCvt cvt;
120 AjPStr subt = NULL;
121
122 ajint b1;
123 ajint b2;
124 ajint e1;
125 ajint e2;
126 AjPStr se1;
127 AjPStr se2;
128 ajint ithresh;
129 AjBool stretch;
130 PPoint ppt = NULL;
131 float xa[1];
132 float ya[1];
133 AjPGraphdata gdata=NULL;
134 AjPStr tit = NULL;
135 AjIList iter = NULL;
136 float x1 = 0.;
137 float x2 = 0.;
138 float y1 = 0.;
139 float y2 = 0.;
140 ajuint tui;
141
142 se1 = ajStrNew();
143 se2 = ajStrNew();
144
145 embInit("dotmatcher", argc, argv);
146
147 seq2 = ajAcdGetSeq("asequence");
148 seq = ajAcdGetSeq("bsequence");
149 stretch = ajAcdGetToggle("stretch");
150 graph = ajAcdGetGraph("graph");
151 xygraph = ajAcdGetGraphxy("xygraph");
152 windowsize = ajAcdGetInt("windowsize");
153 ithresh = ajAcdGetInt("threshold");
154 matrix = ajAcdGetMatrix("matrixfile");
155
156 sub = ajMatrixGetMatrix(matrix);
157 cvt = ajMatrixGetCvt(matrix);
158
159 thresh = (float)ithresh;
160
161 ajtime = ajTimeNew();
162
163 tim = time(0);
164 ajTimeSetLocal(ajtime, tim);
165
166 b1 = ajSeqGetBegin(seq);
167 b2 = ajSeqGetBegin(seq2);
168 e1 = ajSeqGetEnd(seq);
169 e2 = ajSeqGetEnd(seq2);
170 len1 = ajSeqGetLen(seq);
171 len2 = ajSeqGetLen(seq2);
172 tui = ajSeqGetLen(seq);
173 flen1 = (float) tui;
174 tui = ajSeqGetLen(seq2);
175 flen2 = (float) tui;
176
177 ajStrAssignSubC(&se1,ajSeqGetSeqC(seq),b1-1,e1-1);
178 ajStrAssignSubC(&se2,ajSeqGetSeqC(seq2),b2-1,e2-1);
179 ajSeqAssignSeqS(seq,se1);
180 ajSeqAssignSeqS(seq2,se2);
181
182
183 s1 = ajStrGetPtr(ajSeqGetSeqS(seq));
184 s2 = ajStrGetPtr(ajSeqGetSeqS(seq2));
185
186
187 aa0str = ajStrNewRes(1+len1); /* length plus trailing blank */
188 aa1str = ajStrNewRes(1+len2);
189
190 list = ajListNew();
191
192
193 for(i=0;i<len1;i++)
194 ajStrAppendK(&aa0str,(char)ajSeqcvtGetCodeK(cvt, *s1++));
195
196 for(i=0;i<len2;i++)
197 ajStrAppendK(&aa1str,(char)ajSeqcvtGetCodeK(cvt, *s2++));
198
199 max = (float)len1;
200 if(len2 > max)
201 max = (float) len2;
202
203 xmargin = ymargin = max *(float)0.15;
204 ticklen = xmargin*(float)0.1;
205 onefifth = xmargin*(float)0.2;
206
207 subt = ajStrNewC((strret=
208 ajFmtString("(windowsize = %d, threshold = %3.2f %D)",
209 windowsize,thresh,ajtime)));
210
211
212
213 if(!stretch)
214 {
215 if( ajStrGetLen(ajGraphGetSubtitleS(graph)) <=1)
216 ajGraphSetSubtitleS(graph,subt);
217
218 ajGraphOpenWin(graph, (float)0.0-ymargin,(max*(float)1.35)+ymargin,
219 (float)0.0-xmargin,(float)max+xmargin);
220
221 ajGraphicsDrawposTextAtmid(flen1*(float)0.5,
222 (float)0.0-(xmargin/(float)2.0),
223 ajGraphGetXlabelC(graph));
224 ajGraphicsDrawposTextAtlineJustify((float)0.0-(xmargin*(float)0.75),
225 flen2*(float)0.5,
226 (float)0.0-(xmargin*(float)0.75),flen1,
227 ajGraphGetYlabelC(graph),0.5);
228
229 ajGraphicsSetCharscale(0.5);
230 }
231
232
233
234 s1= ajStrGetPtr(aa0str);
235 s2 = ajStrGetPtr(aa1str);
236
237 for(j=0; (ajint)j < (ajint)len2-windowsize;j++)
238 {
239 i =0;
240 total = 0;
241 abovethresh =0;
242
243 k = j;
244 for(l=0;l<windowsize;l++)
245 total = total + sub[(ajint)s1[i++]][(ajint)s2[k++]];
246
247 if(total >= thresh)
248 {
249 abovethresh=1;
250 starti = i-windowsize;
251 startj = k-windowsize;
252 }
253
254 while(i < len1 && k < len2)
255 {
256 total = total - sub[(ajint)s1[i-windowsize]]
257 [(ajint)s2[k-windowsize]];
258 total = total + sub[(ajint)s1[i]][(ajint)s2[k]];
259
260 if(abovethresh)
261 {
262 if(total < thresh)
263 {
264 abovethresh = 0;
265 /* draw the line */
266 dotmatcher_pushpoint(&list,(float)starti,(float)startj,
267 (float)i-1,(float)k-1,stretch);
268 }
269 }
270 else if(total >= thresh)
271 {
272 starti = i-windowsize;
273 startj = k-windowsize;
274 abovethresh= 1;
275 }
276 i++;
277 k++;
278 }
279
280 if(abovethresh)
281 /* draw the line */
282 dotmatcher_pushpoint(&list,(float)starti,(float)startj,
283 (float)i-1,(float)k-1,
284 stretch);
285 }
286
287 for(i=0; (ajint)i < (ajint)len1-windowsize;i++)
288 {
289 j = 0;
290 total = 0;
291 abovethresh =0;
292
293 k = i;
294 for(l=0;l<windowsize;l++)
295 total = total + sub[(ajint)s1[k++]][(ajint)s2[j++]];
296
297 if(total >= thresh)
298 {
299 abovethresh=1;
300 starti = k-windowsize;
301 startj = j-windowsize;
302 }
303
304 while(k < len1 && j < len2)
305 {
306 total = total - sub[(ajint)s1[k-windowsize]]
307 [(ajint)s2[j-windowsize]];
308 total = total + sub[(ajint)s1[k]][(ajint)s2[j]];
309
310 if(abovethresh)
311 {
312 if(total < thresh)
313 {
314 abovethresh = 0;
315 /* draw the line */
316 dotmatcher_pushpoint(&list,(float)starti,(float)startj,
317 (float)k-1,(float)j-1,stretch);
318 }
319 }
320 else if(total >= thresh)
321 {
322 starti = k-windowsize;
323 startj = j-windowsize;
324 abovethresh= 1;
325 }
326 j++;
327 k++;
328 }
329
330 if(abovethresh)
331 /* draw the line */
332 dotmatcher_pushpoint(&list,(float)starti,(float)startj,
333 (float)k-1,(float)j-1,
334 stretch);
335 }
336
337 if(boxit && !stretch)
338 {
339 ajGraphicsDrawposRect(0.0,0.0,flen1, flen2);
340
341 i=0;
342 while(acceptableticksx[i]*numbofticks < len1)
343 i++;
344
345 if(i<=13)
346 tickgap = (float)acceptableticksx[i];
347 else
348 tickgap = (float)acceptableticksx[10];
349 ticklen = xmargin*(float)0.1;
350 onefifth = xmargin*(float)0.2;
351
352 if(len2/len1 > 10 )
353 {
354 /* if a lot smaller then just label start and end */
355 ajGraphicsDrawposLine((float)0.0,(float)0.0,(float)0.0,(float)0.0-ticklen);
356 sprintf(ptr,"%d",b1-1);
357 ajGraphicsDrawposTextAtmid((float)0.0,(float)0.0-(onefifth),ptr);
358
359 ajGraphicsDrawposLine(flen1,(float)0.0,
360 flen1,(float)0.0-ticklen);
361 sprintf(ptr,"%u",len1+b1-1);
362 ajGraphicsDrawposTextAtmid(flen1,(float)0.0-(onefifth),ptr);
363
364 }
365 else
366 for(k2=0.0;k2<len1;k2+=tickgap)
367 {
368 ajGraphicsDrawposLine(k2,(float)0.0,k2,(float)0.0-ticklen);
369 sprintf(ptr,"%d",(ajint)k2+b1-1);
370 ajGraphicsDrawposTextAtmid(k2,(float)0.0-(onefifth),ptr);
371 }
372
373 i = 0;
374 while(acceptableticks[i]*numbofticks < len2)
375 i++;
376
377 tickgap = (float)acceptableticks[i];
378 ticklen = ymargin*(float)0.01;
379 onefifth = ymargin*(float)0.02;
380
381 if(len1/len2 > 10 )
382 {
383 /* if a lot smaller then just label start and end */
384 ajGraphicsDrawposLine((float)0.0,(float)0.0,(float)0.0-ticklen,(float)0.0);
385 sprintf(ptr,"%d",b2-1);
386 ajGraphicsDrawposTextAtend((float)0.0-(onefifth),(float)0.0,ptr);
387
388 ajGraphicsDrawposLine((float)0.0,flen2,(float)0.0-ticklen,
389 flen2);
390 sprintf(ptr,"%u",len2+b2-1);
391 ajGraphicsDrawposTextAtend((float)0.0-(onefifth),flen2,ptr);
392 }
393 else
394 for(k2=0.0;k2<len2;k2+=tickgap)
395 {
396 ajGraphicsDrawposLine((float)0.0,k2,(float)0.0-ticklen,k2);
397 sprintf(ptr,"%d",(ajint)k2+b2-1);
398 ajGraphicsDrawposTextAtend((float)0.0-(onefifth),k2,ptr);
399 }
400 }
401
402
403 if(!stretch)
404 ajGraphicsClose();
405 else /* the xy graph for -stretch */
406 {
407 tit = ajStrNew();
408 ajFmtPrintS(&tit,"%S",ajGraphGetTitleS(xygraph));
409
410
411 gdata = ajGraphdataNewI(1);
412 xa[0] = (float)b1;
413 ya[0] = (float)b2;
414
415 ajGraphSetTitleC(xygraph,ajStrGetPtr(tit));
416
417 ajGraphSetXlabelC(xygraph,ajSeqGetNameC(seq));
418 ajGraphSetYlabelC(xygraph,ajSeqGetNameC(seq2));
419
420 ajGraphdataSetTypeC(gdata,"2D Plot Float");
421 ajGraphdataSetTitleS(gdata,subt);
422 ajGraphdataSetMinmax(gdata,(float)b1,(float)e1,(float)b2,
423 (float)e2);
424 ajGraphdataSetTruescale(gdata,(float)b1,(float)e1,(float)b2,
425 (float)e2);
426 ajGraphxySetXstartF(xygraph,(float)b1);
427 ajGraphxySetXendF(xygraph,(float)e1);
428 ajGraphxySetYstartF(xygraph,(float)b2);
429 ajGraphxySetYendF(xygraph,(float)e2);
430
431 ajGraphxySetXrangeII(xygraph,b1,e1);
432 ajGraphxySetYrangeII(xygraph,b2,e2);
433
434
435 if(list)
436 {
437 iter = ajListIterNewread(list);
438 while((ppt = ajListIterGet(iter)))
439 {
440 x1 = ppt->x1+b1-1;
441 y1 = ppt->y1+b2-1;
442 x2 = ppt->x2+b1-1;
443 y2 = ppt->y2+b2-1;
444 ajGraphAddLine(xygraph,x1,y1,x2,y2,0);
445 AJFREE(ppt);
446 }
447 ajListIterDel(&iter);
448 }
449
450 ajGraphdataAddXY(gdata,xa,ya);
451 ajGraphDataReplace(xygraph,gdata);
452
453
454 ajGraphxyDisplay(xygraph,ajFalse);
455 ajGraphicsClose();
456
457 ajStrDel(&tit);
458 }
459
460
461
462 ajListFree(&list);
463
464 ajSeqDel(&seq);
465 ajSeqDel(&seq2);
466 ajGraphxyDel(&graph);
467 ajGraphxyDel(&xygraph);
468 ajMatrixDel(&matrix);
469 ajTimeDel(&ajtime);
470
471 /* deallocate memory */
472 ajStrDel(&aa0str);
473 ajStrDel(&aa1str);
474 ajStrDel(&se1);
475 ajStrDel(&se2);
476 ajStrDel(&subt);
477
478 AJFREE(strret); /* created withing ajFmtString */
479
480 embExit();
481
482 return 0;
483 }
484
485
486
487
488 /* @funcstatic dotmatcher_pushpoint *******************************************
489 **
490 ** Undocumented.
491 **
492 ** @param [u] l [AjPList*] Undocumented
493 ** @param [r] x1 [float] Undocumented
494 ** @param [r] y1 [float] Undocumented
495 ** @param [r] x2 [float] Undocumented
496 ** @param [r] y2 [float] Undocumented
497 ** @param [r] stretch [AjBool] Do a stretch plot
498 ** @return [void]
499 ** @@
500 ******************************************************************************/
501
dotmatcher_pushpoint(AjPList * l,float x1,float y1,float x2,float y2,AjBool stretch)502 static void dotmatcher_pushpoint(AjPList *l, float x1, float y1, float x2,
503 float y2, AjBool stretch)
504 {
505 PPoint p;
506
507 if(!stretch)
508 {
509 ajGraphicsDrawposLine(x1+1,y1+1,x2+1,y2+1);
510 return;
511 }
512
513 AJNEW0(p);
514 p->x1 = x1+1;
515 p->y1 = y1+1;
516 p->x2 = x2+1;
517 p->y2 = y2+1;
518 ajListPush(*l,(void *)p);
519
520 return;
521 }
522