1 /******************************************************************************
2
3 #### ##### #### #### #### ## # # ####
4 # # # # # # # # # # ## # # #
5 # # # # #### #### # # # # # # #
6 # # # ##### # # # ###### # # # ### #
7 # # # # # # # # # # # # ## ### # #
8 ### # # #### ####### #### #### # # # # ### ####
9
10 ******************************************************************************/
11 /* This file is part of MAPMAKER 3.0b, Copyright 1987-1992, Whitehead Institute
12 for Biomedical Research. All rights reserved. See READ.ME for license. */
13
14 #define INC_LIB
15 #define INC_SHELL
16 #define INC_CALLQCTM
17 #define INC_QTOPLEVEL
18 #define INC_QLOWLEVEL
19 #include "qtl.h"
20
21 void ps_file_start(), ps_file_end(), ps_page_start(), ps_page_end();
22 void do_bezier(), draw_axes(), draw_x();
23
24 #define XZERO 0.0
25 #define YZERO 0.0
26
27 #define SOLID_LINE ""
28 #define THICK_LINE "2 LW"
29 #define THIN_LINE ".5 LW"
30 #define EVEN_DASH "[4 4] 0 setdash"
31 #define BIG_DASH "[6 2] 0 setdash"
32 #define SMALL_DASH "[2 6] 0 setdash"
33 #define DOTTED "[1 2] 0 setdash"
34
35 real xscale, yscale, best;
36 #define MAX_CHROM_LEN 1000
37
print_ps_wiggle_order(wiggle,order,threshold)38 void print_ps_wiggle_order(wiggle, order, threshold)
39 int wiggle, order;
40 real threshold;
41 {
42 int i, j, k, count=0, num_notches=0, pagenum=1;
43 double *xval, *yval, *notch, highest, longest, current_len;
44 WIGGLE_OPERATION *op;
45 WIGGLE_INTERVAL **data;
46 WIGGLE_POINT *point;
47 FILE *fp;
48 char **label, *trait_str, *filename;
49
50 best = 0.0;
51 array(xval, MAX_CHROM_LEN, real);
52 array(yval, MAX_CHROM_LEN, real);
53 array(notch, MAX_CHROM_LOC, real);
54 matrix(label, MAX_CHROM_LOC, NAME_LEN+1, char);
55 filename = get_temp_string();
56 sf(filename,"scan%d_%d.ps",wiggle+1,order+1);
57 fp = fopen(filename,"w");
58
59 ps_file_start(fp);
60 ps_page_start(fp, pagenum);
61
62 op = wiggles[wiggle]; data = op->data[order];
63
64 /***** Read through the data and calculate scale based on *****/
65 /***** longest chromosome and highest lodscore. *****/
66
67 longest=0.0; highest=0.0; current_len=0.0;
68 for (i=0; i<op->num_wiggled_intervals; i++) {
69 if (i>0 && !data[i]->contig) {
70 /* chromosome done */
71 if(current_len > longest) longest = current_len;
72 current_len = 0.0;
73 }
74 for (j=0; j<data[i]->num_points; j++) {
75 point=data[i]->point[j];
76 if(point->lod_score > highest) highest = point->lod_score;
77 current_len += 2.0;
78 }
79 }
80 if(current_len > longest) longest = current_len;
81 xscale = (double) 700.0/longest;
82 yscale = rminf(50.0, 400.0/highest);
83
84 trait_str = get_temp_string();
85 sf(trait_str,"LOD score - Trait %d (%s)", op->trait+1,
86 raw.trait_name[op->trait]);
87 if(!print_names) sf(label[0],"%d",
88 raw.original_locus[data[0]->map->left[data[0]->map->num_intervals-1]]);
89 else sf(label[0],"%s",
90 raw.locus_name[data[0]->map->left[data[0]->map->num_intervals-1]]);
91 for (i=0; i<op->num_wiggled_intervals; i++) {
92 if (i>0) {
93 notch[num_notches++] = count*2.0;
94 if (!print_names) sf(label[num_notches],"%d",
95 raw.original_locus[data[i-1]->map->right[data[i-1]->map->num_intervals-1]]);
96 else sf(label[num_notches],"%s",
97 raw.locus_name[data[i-1]->map->right[data[i-1]->map->num_intervals-1]]);
98 }
99 if (i>0 && !data[i]->contig) {
100 /* a page is done */
101 draw_axes(fp, notch, num_notches, label, trait_str, threshold);
102 do_bezier(fp, xval, yval, count, 0.0, 0.0, SOLID_LINE);
103 ps_page_end(fp);
104 pagenum++;
105 ps_page_start(fp, pagenum);
106 if(!print_names) sf(label[0],"%d",
107 raw.original_locus[data[i]->map->left[data[i]->map->num_intervals-1]]);
108 else sf(label[0],"%s",
109 raw.locus_name[data[i]->map->left[data[i]->map->num_intervals-1]]);
110 count = 0; num_notches = 0;
111 best = 0.0;
112 }
113 for(j=0; j<data[i]->num_points; j++) {
114 point=data[i]->point[j];
115 xval[count] = count*2.0;
116 yval[count] = point->lod_score;
117 if(point->lod_score>best) best = point->lod_score;
118 count++;
119 }
120 }
121
122 /* dump last page */
123 notch[num_notches++] = count*2.0;
124 if(!print_names) sf(label[num_notches],"%d",
125 raw.original_locus[data[i-1]->map->right[data[i-1]->map->num_intervals-1]]);
126 else sf(label[num_notches],"%s",
127 raw.locus_name[data[i-1]->map->right[data[i-1]->map->num_intervals-1]]);
128
129 draw_axes(fp, notch, num_notches, label, trait_str, threshold);
130 do_bezier(fp, xval, yval, count, 0.0, 0.0, SOLID_LINE);
131
132 ps_page_end(fp);
133 ps_file_end(fp);
134
135 fclose(fp);
136 sf(ps,"scan %d.%d saved in PostScript file '%s'\n",
137 wiggle+1,order+1,filename); pr();
138 unarray(xval, real);
139 unarray(yval, real);
140 unarray(notch, real);
141 unmatrix(label, MAX_CHROM_LOC, char);
142 }
143
144
draw_axes(fp,xnotch,num_notches,label,y_name,dotted_val)145 void draw_axes(fp, xnotch, num_notches, label, y_name, dotted_val)
146 FILE *fp;
147 double *xnotch;
148 int num_notches;
149 char **label, *y_name;
150 double dotted_val;
151 {
152 int i;
153 double prev, next, current;
154
155 fprintf(fp,"500 50 translate\n");
156 fprintf(fp,"90 rotate\n");
157 fprintf(fp,"%.2lf %.2lf moveto\n", XZERO, YZERO);
158 fprintf(fp,"GS 0.5 LW 700 0 rlineto stroke GR\n");
159 fprintf(fp,"GS 0.5 LW 0 450 rlineto stroke GR\n");
160 fprintf(fp,"GS 0 460 rmoveto /Times-Roman FF 14 SF F (%s)S GR\n",y_name);
161
162 for(i=1; i<=(int)(best+1.0); i++) {
163 fprintf(fp,"%.2lf %.2lf moveto\n", XZERO-4, YZERO+((double)i*yscale));
164 fprintf(fp,"GS 0.5 LW 9 0 rlineto stroke GR\n");
165 fprintf(fp,"GS -3 -3 rmoveto /Times-Roman FF 9 SF F (%d.0)SR GR\n",i);
166 }
167 /* draw dotted line - y = dotted_val */
168 fprintf(fp,"%.2lf %.2lf moveto\n", XZERO, YZERO+(dotted_val*yscale));
169 fprintf(fp,"GS 0.5 LW [4 4] 0 setdash 700 0 rlineto stroke GR\n");
170
171 prev = 0.0;
172 for (i=0; i<num_notches; i++) {
173 current = XZERO+(xnotch[i]*xscale);
174 next = (i < num_notches-1) ? XZERO+(xnotch[i+1]*xscale) : 999.9;
175
176 if(next-current < 5.0) {
177 if(current-prev > 6.0) current -= 1.0;
178 }
179 fprintf(fp,"%.2lf %.2lf moveto\n",current,YZERO-4);
180 fprintf(fp,"GS 0.5 LW 0 2 rmoveto 0 4 rlineto stroke GR\n");
181 fprintf(fp,"GS -90 rotate 0 -3 rmoveto /Times-Roman FF 8 SF F (%s)S GR\n",
182 label[i+1]);
183 prev = current;
184 }
185
186 fprintf(fp,"%.2lf %.2lf moveto\n", XZERO, YZERO-4);
187 fprintf(fp,"GS -90 rotate 0 -3 rmoveto /Times-Roman FF 8 SF F (%s)S GR\n",
188 label[0]);
189 draw_x(fp);
190 }
191
192
do_bezier(fp,xval,yval,num_points,s0,sn,line_type)193 void do_bezier(fp, xval, yval, num_points, s0, sn, line_type)
194 FILE *fp;
195 double *xval, *yval, s0, sn;
196 int num_points;
197 char *line_type;
198 {
199 int i, j;
200 double *x0, *x1, *x2, *x3, *y0, *y1, *y2, *y3, *slope, s1, s2;
201
202 array(x0, num_points, real);
203 array(x1, num_points, real);
204 array(x2, num_points, real);
205 array(x3, num_points, real);
206 array(y0, num_points, real);
207 array(y1, num_points, real);
208 array(y2, num_points, real);
209 array(y3, num_points, real);
210 array(slope, num_points, real);
211
212 /* Bezier functionality draws an appropriate curve between points (x0,y0)
213 and (x3,y3) based on Bezier cubic control points (x1,y1) and (x2,y2)
214 (x1,y1) lies on the tangent from (x0,y0) and (x2,y2) lies on the tangent
215 from (x3,y3) with x0 < x1 < x2 < x3 for our cases since
216 x3 = x0 + scan_interval_length (usually 2 cM) */
217
218
219 /* (x0,y0) and (x3,y3) are the known data points, scaled to size of page */
220
221 for(i = 0; i < num_points-1; i++) {
222 x0[i] = xval[i]*xscale + XZERO;
223 y0[i] = yval[i]*yscale + YZERO;
224 x3[i] = xval[i+1]*xscale + XZERO;
225 y3[i] = yval[i+1]*yscale + YZERO;
226 }
227
228 /* estimate the slope of the tangent at each data point...given three
229 points on a curve (x1,y1) (x2,y2) (x3,y3), a good guess for the
230 slope of the tangent at (x2,y2) is the average of the slopes of the
231 lines connecting (x1,y1)(x2,y2) and (x2,y2)(x3,y3) - since this method
232 does not work for the first or last point on the curve for obvious
233 reasons, those values are handed into the function as s0 and sn...
234 usually values of 0.0 will be just fine */
235
236 slope[0] = s0;
237 for(i = 1; i < num_points-1; i++) {
238 s1 = (y3[i-1]-y0[i-1])/(x3[i-1]-x0[i-1]);
239 s2 = (y3[i]-y0[i])/(x3[i]-x0[i]);
240 slope[i] = (s1+s2)/2.0;
241 if(slope[i] > 10.0) slope[i] = 10.0;
242 if(slope[i] < -10.0) slope[i] = -10.0;
243 }
244 slope[i] = sn;
245
246
247 /* calculate points (x1,y1) and (x2,y2) now that we have approximate
248 tangent lines for our real data points */
249
250 for(i = 0; i < num_points-1; i++) {
251 x1[i] = (x3[i]-x0[i])/3.0;
252 y1[i] = y0[i] + x1[i]*slope[i];
253 x2[i] = x1[i];
254 y2[i] = y3[i] - x2[i]*slope[i+1];
255 x1[i] += x0[i];
256 x2[i] += x1[i];
257 }
258
259 /* Postscript conveniently draws a Bezier cubic section from the current
260 position to position (x3,y3) based on the calculated control points
261 (x1,y1) and (x2,y2) */
262
263 fprintf(fp,"%.2lf %.2lf moveto\n",x0[0],y0[0]);
264
265 for(i = 0; i < num_points-1; i++) {
266 fprintf(fp,"GS %s %.2lf %.2lf %.2lf %.2lf %.2lf %.2lf curveto stroke GR\n",
267 line_type, x1[i], y1[i], x2[i], y2[i], x3[i], y3[i]);
268 fprintf(fp,"%.2lf %.2lf moveto\n", x3[i], y3[i]);
269 draw_x(fp);
270 }
271 fprintf(fp,"%.2lf %.2lf moveto\n", x3[i-1], y3[i-1]);
272 draw_x(fp);
273
274 unarray(x0, real);
275 unarray(x1, real);
276 unarray(x2, real);
277 unarray(x3, real);
278 unarray(y0, real);
279 unarray(y1, real);
280 unarray(y2, real);
281 unarray(y3, real);
282 unarray(slope, real);
283 }
284
285
draw_x(fp)286 void draw_x(fp)
287 FILE *fp;
288 {
289 #ifdef DRAWX
290 fprintf(fp,"GS -3 -3 rlineto stroke GR\n");
291 fprintf(fp,"GS 3 3 rlineto stroke GR\n");
292 fprintf(fp,"GS 3 -3 rlineto stroke GR\n");
293 fprintf(fp,"GS -3 3 rlineto stroke GR\n");
294 #endif
295 }
296
297
298
ps_file_start(fp)299 void ps_file_start(fp)
300 FILE *fp;
301 {
302 fprintf(fp,"%%!PS-Adobe-3.0\n");
303 fprintf(fp,"%%%%Creator: MAPMAKER\n");
304 fprintf(fp,"%%%%LanguageLevel: 1\n");
305 fprintf(fp,"%%%%PageOrder: Special\n");
306
307 fprintf(fp,"%%%%BeginProlog\n");
308 fprintf(fp,"%%%%BeginResource: procset Map_Painter_prolog\n");
309 fprintf(fp,"/Map_Painter_prolog 100 dict def\n");
310 fprintf(fp,"Map_Painter_prolog begin\n");
311 fprintf(fp,"/CF {dup 0 eq\n");
312 fprintf(fp," {pop /Times-Bold}\n");
313 fprintf(fp," {1 eq {/Times-Roman} {/Times-Italic} ifelse}\n");
314 fprintf(fp," ifelse} def\n");
315 fprintf(fp,"/F {setfont} def\n");
316 fprintf(fp,"/FF {findfont} def\n");
317 fprintf(fp,"/GM {restore} def\n");
318 fprintf(fp,"/GR {grestore} def\n");
319 fprintf(fp,"/GS {gsave} def\n");
320 fprintf(fp,"/LC {currentpoint pop dup /l exch def add /r exch def\n");
321 fprintf(fp," counttomark 3 idiv {\n");
322 fprintf(fp," currentpoint pop r ge {l currentpoint exch pop LLD sub moveto} if\n");
323 fprintf(fp," LD\n");
324 fprintf(fp," } repeat pop} def\n");
325 fprintf(fp,"/LD {/f exch def /n exch def /p exch def\n");
326 fprintf(fp," p length 0 ne {\n");
327 fprintf(fp," /Times-Roman findfont LFS 1 sub scalefont setfont p show\n");
328 fprintf(fp," 2.5 0 rmoveto\n");
329 fprintf(fp," } if\n");
330 fprintf(fp," f CF findfont LFS scalefont setfont\n");
331 fprintf(fp," n show 6 0 rmoveto} def\n");
332 fprintf(fp,"/LM {save} def\n");
333 fprintf(fp,"/LW {setlinewidth} def\n");
334 fprintf(fp,"/S {show} def\n");
335 fprintf(fp,"/SF {scalefont} def\n");
336 fprintf(fp,"/SR {dup stringwidth pop neg 0 rmoveto show} def\n");
337 fprintf(fp,"/TR {translate} def\n");
338 fprintf(fp,"/XY {moveto} def\n");
339 fprintf(fp,"/VSdict 4 dict def\n");
340 fprintf(fp,"/VS { VSdict begin\n");
341 fprintf(fp," /thestring exch def\n");
342 fprintf(fp," /lineskip exch def\n");
343 fprintf(fp," thestring\n");
344 fprintf(fp," {\n");
345 fprintf(fp," /charcode exch def\n");
346 fprintf(fp," /thechar ( ) dup 0 charcode put def\n");
347 fprintf(fp," 0 lineskip neg rmoveto\n");
348 fprintf(fp," gsave\n");
349 fprintf(fp," thechar stringwidth pop 2 div neg 0 rmoveto\n");
350 fprintf(fp," thechar show\n");
351 fprintf(fp," grestore\n");
352 fprintf(fp," } forall\n");
353 fprintf(fp," end\n");
354 fprintf(fp," } def\n");
355
356 fprintf(fp,"end\n");
357 fprintf(fp,"%%%%EndResource\n");
358 fprintf(fp,"%%%%EndProlog\n");
359
360 fprintf(fp,"/LFS 8 def /LLD 7.5 def\n");
361 fprintf(fp,"%%%%BeginSetup\n");
362 fprintf(fp,"Map_Painter_prolog begin\n");
363 fprintf(fp,"gsave\n");
364 fprintf(fp,"%%%%EndSetup\n");
365 }
366
ps_file_end(fp)367 void ps_file_end(fp)
368 FILE *fp;
369 {
370 fprintf(fp,"%%%%Trailer\n");
371 fprintf(fp,"grestore\n");
372 fprintf(fp,"end %% Map_Painter_prolog\n");
373 fprintf(fp,"%%%%EOF\n");
374 }
375
ps_page_start(fp,pagenum)376 void ps_page_start(fp,pagenum)
377 FILE *fp;
378 int pagenum;
379 {
380 fprintf(fp,"%%%%Page: ? %d\n",pagenum);
381 fprintf(fp,"%%%%BeginPageSetup\n");
382 fprintf(fp,"LM\n");
383 fprintf(fp,"%%%%EndPageSetup\n");
384 }
385
ps_page_end(fp)386 void ps_page_end(fp)
387 FILE *fp;
388 {
389 fprintf(fp,"GM showpage\n");
390 }
391
392
393
394
print_ps_multi_wiggle(wiggle,threshold)395 void print_ps_multi_wiggle(wiggle, threshold)
396 int wiggle;
397 real threshold;
398 {
399 int i, j, k, count=0, num_notches=0, pagenum=1, order;
400 double **xval, **yval, *notch, highest, longest, current_len;
401 WIGGLE_OPERATION *op;
402 WIGGLE_INTERVAL **data;
403 WIGGLE_POINT *point;
404 FILE *fp;
405 char **label, *trait_str, *filename, *line_choice();
406
407 op = wiggles[wiggle];
408
409 if(op->num_orders > 4) error("too many wiggle orders to display on one graph");
410
411 best = 0.0;
412 matrix(xval, op->num_orders, MAX_CHROM_LEN, real);
413 matrix(yval, op->num_orders, MAX_CHROM_LEN, real);
414 array(notch, MAX_CHROM_LOC, real);
415 matrix(label, MAX_CHROM_LOC, NAME_LEN+1, char);
416 filename = get_temp_string();
417 sf(filename,"scan%d_x.ps",wiggle+1);
418 fp = fopen(filename,"w");
419
420 ps_file_start(fp);
421 ps_page_start(fp, pagenum);
422 longest=0.0; highest=0.0; current_len=0.0;
423 for (i=0; i<op->num_wiggled_intervals; i++) {
424 if (i>0 && !op->data[0][i]->contig) {
425 if(current_len > longest) longest = current_len;
426 current_len = 0.0;
427 }
428 for (j=0; j<op->data[0][i]->num_points; j++) {
429 point=op->data[0][i]->point[j];
430 if(point->lod_score > highest) highest = point->lod_score;
431 current_len += 2.0;
432 }
433 }
434
435 if(current_len > longest) longest = current_len;
436 xscale = (double) 700.0/longest;
437 yscale = rminf(50.0, 400.0/highest);
438
439 trait_str = get_temp_string();
440 sf(trait_str,"LOD score - Trait %d (%s)", op->trait+1,
441 raw.trait_name[op->trait]);
442 if (!print_names) sf(label[0],"%d",
443 raw.original_locus[op->data[0][0]->map->left[op->data[0][0]->map->num_intervals-1]]);
444 else sf(label[0],"%s",
445 raw.locus_name[op->data[0][0]->map->left[op->data[0][0]->map->num_intervals-1]]);
446 for (i=0; i<op->num_wiggled_intervals; i++) {
447 if (i>0) {
448 notch[num_notches++] = count*2.0;
449 if (!print_names) sf(label[num_notches],"%d",
450 raw.original_locus[op->data[0][i-1]->map->right[op->data[0][i-1]->map->num_intervals-1]]);
451 else sf(label[num_notches],"%s",
452 raw.locus_name[op->data[0][i-1]->map->right[op->data[0][i-1]->map->num_intervals-1]]);
453 }
454 if (i>0 && !op->data[0][i]->contig) {
455
456 draw_axes(fp, notch, num_notches, label, trait_str, threshold);
457 for (order=0; order<op->num_orders; order++)
458 do_bezier(fp, xval[order], yval[order], count, 0.0, 0.0, line_choice(order));
459
460 fprintf(fp,"0 -40 moveto\n");
461 fprintf(fp,"GS /Times-Roman FF 9 SF F (GENETICS:)S GR\n");
462 for (order=0; order<op->num_orders; order++) {
463 fprintf(fp,"0 %d moveto\n",-50-(order*10));
464 fprintf(fp,"GS /Times-Roman FF 9 SF F (%s)S GR\n",
465 genetics_str(&op->data[order][0]->map->constraint[op->num_intervals-1],FALSE));
466 fprintf(fp,"%d %d moveto\n",50,-50-(order*10));
467 fprintf(fp,"GS %s 20 0 rlineto stroke GR\n",line_choice(order));
468 }
469 ps_page_end(fp);
470 pagenum++;
471 ps_page_start(fp, pagenum);
472 if (!print_names) sf(label[0],"%d",
473 raw.original_locus[op->data[0][i]->map->left[op->data[0][i]->map->num_intervals-1]]);
474 else sf(label[0],"%s",
475 raw.locus_name[op->data[0][i]->map->left[op->data[0][i]->map->num_intervals-1]]);
476 count = 0; num_notches = 0;
477 best = 0.0;
478 }
479 for (j=0; j<op->data[0][i]->num_points; j++) {
480 for (order=0; order<op->num_orders; order++) {
481 point=op->data[order][i]->point[j];
482 xval[order][count] = count*2.0;
483 yval[order][count] = point->lod_score;
484 if(point->lod_score>best) best = point->lod_score;
485 }
486 count++;
487 }
488 }
489
490 notch[num_notches++] = count*2.0;
491 if (!print_names) sf(label[num_notches],"%d",
492 raw.original_locus[op->data[0][i-1]->map->right[op->data[0][i-1]->map->num_intervals-1]]);
493 else sf(label[num_notches],"%s",
494 raw.locus_name[op->data[0][i-1]->map->right[op->data[0][i-1]->map->num_intervals-1]]);
495 draw_axes(fp, notch, num_notches, label, trait_str, threshold);
496 for (order=0; order<op->num_orders; order++)
497 do_bezier(fp, xval[order], yval[order], count, 0.0, 0.0,
498 line_choice(order));
499
500 fprintf(fp,"0 -40 moveto\n");
501 fprintf(fp,"GS /Times-Roman FF 9 SF F (GENETICS:)S GR\n");
502 for (order=0; order<op->num_orders; order++) {
503 fprintf(fp,"0 %d moveto\n",-50-(order*10));
504 fprintf(fp,"GS /Times-Roman FF 9 SF F (%s)S GR\n",
505 genetics_str(&op->data[order][0]->map->constraint[op->num_intervals-1],
506 FALSE));
507 fprintf(fp,"%d %d moveto\n",50,-50-(order*10));
508 fprintf(fp,"GS %s 20 0 rlineto stroke GR\n",line_choice(order));
509 }
510
511 ps_page_end(fp);
512 ps_file_end(fp);
513
514 fclose(fp);
515 sf(ps,"scan %d.x saved in PostScript file '%s'\n",wiggle+1,filename); pr();
516 unarray(xval, real);
517 unarray(yval, real);
518 unarray(notch, real);
519 unmatrix(label, MAX_CHROM_LOC, char);
520 }
521
line_choice(order)522 char *line_choice(order)
523 int order;
524 {
525 switch(order) {
526 case 0: return(THICK_LINE);
527 case 1: return(THIN_LINE);
528 case 2: return(DOTTED);
529 case 3: return(EVEN_DASH);
530 default: return(SOLID_LINE);
531 }
532 }
533
534