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 /********** QPRINT.C - QTL's PRINTING ROUTINES **********/
15
16 #define INC_LIB
17 #define INC_SHELL
18 #define INC_CALLQCTM
19 #define INC_QTOPLEVEL
20 #define INC_QLOWLEVEL /* why? */
21 #include "qtl.h"
22
23 char *interval_str();
24 char *trait_str();
25
26 void print_test_wiggle_map();
27 void map_printer();
28
29
30 /* QTL MAP OUTPUT FORMAT: STEVE changed! Again!
31
32 ==============================================================
33 QTL map for trait 1 (no_error):
34
35 INTERVALS LENGTH QTL-POS WEIGHT DOMINANCE
36 234-678 100.0 100.0 -0.000 -0.000
37 1-2 100.0 100.0 0.000 0.000
38
39 CONTINUOUS-VARIABLE WEIGHT
40 123 (1234567890) -0.000
41
42 INTERVALS LENGTH POSITION GENETICS WEIGHT DOMINANCE <- !print_names
43 234-678 100.0 100.0 recessive -0.000 -0.000
44 1-2 100.0 100.0 free 0.000 0.000
45
46 INTERVALS LENGTH POSITION GENETICS WEIGHT DOMINANCE
47 1234567890 1234567890 100.0 100.0 recessive 0.000 0.000
48
49 chi^2= 12.456 (12 D.F.) log-likelihood= 1234.67
50 mean= 12.456 sigma^2= 12.456 variance-explained= 0.1 %
51 ==============================================================
52 */
53
54 #define INT_TITLE_LINE "INTERVAL%s%s LENGTH QTL-POS "
55 #define INT_LINE "%s %-5s %-5s "
56 #define GENETICS_TITLE "GENETICS "
57 #define GENETICS_FORMAT " %-9s "
58 #define WEIGHT_TITLE "WEIGHT "
59 #define WEIGHT_FORMAT "%s "
60 #define DOM_TITLE "DOMINANCE"
61 #define DOM_FORMAT "%s"
62 #define LOD_TITLE " %VAR LOD" /* FROB */
63 #define LOD_FORMAT " %4.2lf%% %-5.2lf"
64 #define CV_TITLE "CONTINUOUS-VARIABLE WEIGHT\n"
65 #define CV_LINE "%-20s %s\n"
66
67 #define LIKE_LINE \
68 "chi^2= %-6.3lf (%-2d D.F.) log-likelihood= %-7.2lf\n"
69 #define VAR_LINE \
70 "mean= %-6.3lf sigma^2= %-6.3lf variance-explained= %-4.1lf%%\n"
71
72 #define ITER_LINE \
73 "log-likelihood= %-10.5lf delta-log-like= %-10.5lf\n"
74 #define EXTRA_LINE1 \
75 "abs-log-like= %-10.5lf no-data-like= %-10.5lf\n"
76 #define EXTRA_LINE2 \
77 "null-log-like= %-10.5lf like-diff= %-10.5lf\n"
78
79 #define SHORT_MAP_MAX_STARS 10
80
print_map_divider()81 void print_map_divider()
82 {
83 if (print_names && raw.data_type==INTERCROSS) print(LONG_MAP_DIVIDER);
84 else print(MAP_DIVIDER);
85 }
86
87
print_qtl_map(map,free_genetics)88 void print_qtl_map(map,free_genetics) /* has CONT_VAR hooks */
89 QTL_MAP *map;
90 bool *free_genetics;
91 {
92 int i, print_genetics, df_per;
93
94 if (free_genetics==NULL) print_genetics=TRUE;
95 else for (i=0, print_genetics=FALSE; i<map->num_intervals; i++)
96 if (!free_genetics[i]) { print_genetics=TRUE; break; }
97
98 /* hold(6+map->num_intervals); */
99 map_printer(map,print_genetics); nl();
100 df_per= (raw.data_type==INTERCROSS ? 2:1);
101 sf(ps,LIKE_LINE,map->chi_sq,df_per*map->num_intervals,map->log_like); pr();
102 sf(ps,VAR_LINE,map->mu,map->sigma_sq,map->var_explained*100.0); pr();
103
104 /* sf(ps,EXTRA_LINE1,map->abs_log_like,map->no_data_like); pr();
105 sf(ps,EXTRA_LINE2,map->null_log_like,
106 map->no_data_like-map->null_log_like); pr(); */
107 }
108
109
map_printer(map,print_genetics)110 void map_printer(map,print_genetics)
111 QTL_MAP *map;
112 bool print_genetics;
113 {
114 int i;
115 char *trait_string= get_temp_string();
116
117 if (map->num_intervals>0) {
118 sf(ps,INT_TITLE_LINE,map->num_intervals>1 ? "S":" ",
119 print_names ? " ":""); pr();
120 if (print_genetics) print(GENETICS_TITLE);
121 print(WEIGHT_TITLE);
122 if (raw.data_type==INTERCROSS) print(DOM_TITLE);
123 nl();
124 }
125 for (i=0; i<map->num_intervals; i++) {
126 sf(ps,INT_LINE,interval_str(map->left[i],map->right[i],TRUE),
127 dist_str(map->interval_len[i],FALSE),
128 dist_str(map->qtl_pos[i],FALSE)); pr();
129 if (print_genetics) {
130 sf(ps,GENETICS_FORMAT,genetics_str(&map->constraint[i],FALSE));
131 pr();
132 }
133 sf(ps,WEIGHT_FORMAT,rsn(7.5,map->qtl_weight[i])); pr();
134 if (raw.data_type==INTERCROSS)
135 { sf(ps,DOM_FORMAT,rsn(7.5,map->qtl_dominance[i])); pr(); }
136 nl();
137 }
138
139 if (map->num_continuous_vars>0) {
140 if (map->num_intervals>0) nl(); /* divider line */
141 print(CV_TITLE);
142 }
143 for (i=0; i<map->num_continuous_vars; i++) {
144 if (map->cont_var[i]==EPISTASIS_TERM) strcpy(trait_string,"epistasis");
145 else sf(trait_string,"%d (%s)",map->cont_var[i]+1,
146 raw.trait_name[map->cont_var[i]]);
147
148 sf(ps,CV_LINE,trait_string,rsn(7.5,map->cont_var_weight[i])); pr();
149 }
150 }
151
152
print_short_title()153 void print_short_title() /* FROB */
154 {
155 sf(ps,INT_TITLE_LINE," ",print_names ? " ":""); pr();
156 print(WEIGHT_TITLE);
157 if (raw.data_type==INTERCROSS) print(DOM_TITLE);
158 print(LOD_TITLE);
159 nl();
160 }
161
162
print_short_qtl_map(map,threshold,scale)163 void print_short_qtl_map(map,threshold,scale) /* FROB */
164 QTL_MAP *map;
165 real threshold, scale;
166 {
167 int i, n, last;
168 last= map->num_intervals-1;
169
170 sf(ps,INT_LINE,interval_str(map->left[last],map->right[last],TRUE),
171 dist_str(map->interval_len[last],FALSE),
172 dist_str(map->qtl_pos[last],FALSE)); pr();
173 sf(ps,WEIGHT_FORMAT,rsn(7.5,map->qtl_weight[last])); pr();
174 if (raw.data_type==INTERCROSS)
175 { sf(ps,DOM_FORMAT,rsn(7.5,map->qtl_dominance[last])); pr(); }
176 sf(ps,LOD_FORMAT,map->var_explained*100.0,map->log_like); pr();
177 if (map->log_like>=threshold) {
178 n= ((int)((map->log_like-threshold)/scale))+1;
179 if (n>SHORT_MAP_MAX_STARS) n=SHORT_MAP_MAX_STARS;
180 for (i=0; i<n; i++) print("*");
181 }
182 nl();
183 }
184
185
print_iteration(iter,map,delta_log_like)186 void print_iteration(iter,map,delta_log_like)
187 int iter;
188 QTL_MAP *map;
189 real delta_log_like;
190 {
191 print(ITER_DIVIDER);
192 sf(ps,"iteration #%d:\n",iter+1); print(ps);
193 map_printer(map,FALSE);
194 sf(ps,ITER_LINE,map->log_like,delta_log_like);
195 print(ps);
196 sf(ps,VAR_LINE,map->mu,map->sigma_sq,map->var_explained*100.0);
197 print(ps);
198 }
199
200
201 #define NULL_MAP_LINE \
202 "null-like= %-8.3lf null_sigma^2= %-6.3lf \nno-data-like=%-8.3lf\n"
203
print_null_iteration(map)204 void print_null_iteration(map)
205 QTL_MAP *map;
206 {
207 print(ITER_DIVIDER);
208 sf(ps,NULL_MAP_LINE,map->null_log_like,map->null_sigma_sq,
209 map->no_data_like);
210 print(ps);
211 }
212
213
214 /* has CONT_VAR hooks */
do_print_E_step(expected_genotype,S_matrix,expected_recs,n_individuals,n_genotype_vars,n_continuous_vars,n_intervals)215 void do_print_E_step(expected_genotype,S_matrix,expected_recs,n_individuals,
216 n_genotype_vars,n_continuous_vars,n_intervals)
217 real **expected_genotype, **S_matrix;
218 GENO_PROBS *expected_recs;
219 int n_individuals, n_genotype_vars, n_continuous_vars, n_intervals;
220 {
221 int i, j, n, g, q, m, L, R, num;
222
223 num= n_genotype_vars + n_continuous_vars;
224
225 nl();
226 print("expected_genotype:\nindiv");
227 for (n=0; n<num+1; n++)
228 { sf(ps," %-2d ",n); pr(); } nl();
229 for (i=0; i<n_individuals; i++) {
230 sf(ps,"%-4d ",i); pr();
231 for (n=0; n<num+1; n++)
232 { sf(ps,"%-6.3lf ",expected_genotype[i][n]); pr(); }
233 nl();
234 }
235
236 print("\nS_matrix:\n ");
237 for (n=0; n<num+1; n++) { sf(ps," %-2d ",n); pr(); } nl();
238 for (n=0; n<num+1; n++) {
239 sf(ps,"%-2d ",n); pr();
240 for (m=0; m<num+1; m++)
241 { sf(ps,"%-6.2lf ",S_matrix[n][m]); pr(); }
242 nl();
243 }
244
245 print("\nexpected recs: ");
246 for (j=0; j<n_intervals; j++) {
247 sf(ps,"interval %d:\n",j); pr();
248 for_locus_genotypes(raw.data_type,q) {
249 for_interval_genotypes(raw.data_type,g)
250 { sf(ps,"[%d][%d] ",g,q); pr(); } nl();
251 for_interval_genotypes(raw.data_type,g)
252 { sf(ps,"%7.5lf ",expected_recs[j][g][q]); pr(); } nl();
253 }
254 }
255 nl();
256 }
257
258
259 /* WIGGLES OUTPUT FORMAT.... (B1 AND F2, RESPECTIVELY)
260
261 POS WEIGHT %VAR LOG-LIKE |
262 -------------------------------| (1 2) 10.0 cm
263 123.5 -2.456 12.4% 123.567 | 1234567890123456789012345678901234567890123456
264 -------------------------------|
265
266 POS WEIGHT DOM %VAR LOG-LIKE |
267 ---------------------------------------| (1234567890 1234567890) 100.0 cm
268 123.5 -2.456 -2.456 12.4% 123.567 | 12345678901234567890123456789012345678
269 ---------------------------------------|
270 */
271
272 #define B1_WIGGLE_TITLE "POS WEIGHT %VAR LOG-LIKE | "
273 #define B1_WIGGLE_MAP "%5s %s %4.1lf%% %s | "
274 #define B1_NULL_WIGGLE "-------------------------------|\n"
275 #define B1_WIGGLE_INTERVAL \
276 "-------------------------------| %s %s %s\n"
277 #define B1_SPACES_LEFT 43
278
279 #define F2_WIGGLE_TITLE "POS WEIGHT DOM %VAR LOG-LIKE | "
280 #define F2_WIGGLE_MAP "%5s %s %s %4.1lf%% %s | "
281 #define F2_NULL_WIGGLE "---------------------------------------|\n"
282 #define F2_WIGGLE_INTERVAL \
283 "---------------------------------------| %s %s %s\n"
284 #define F2_SPACES_LEFT 35
285
286
print_wiggle_title()287 void print_wiggle_title()
288 {
289 if (raw.data_type==BACKCROSS) { print(B1_WIGGLE_TITLE); }
290 else { print(F2_WIGGLE_TITLE); }
291 nl();
292 }
293
294
print_wiggle_interval(map)295 void print_wiggle_interval(map)
296 QTL_MAP *map;
297 {
298 int i;
299
300 if (map==NULL) {
301 if (raw.data_type==BACKCROSS) strcpy(ps,B1_NULL_WIGGLE);
302 else strcpy(ps,F2_NULL_WIGGLE);
303 } else {
304 i= map->num_intervals-1;
305 sf(ps,
306 (raw.data_type==BACKCROSS ? B1_WIGGLE_INTERVAL:F2_WIGGLE_INTERVAL),
307 interval_str(map->left[i],map->right[i],FALSE),
308 dist_str(map_length(map->left[i],map->right[i]),FALSE),
309 units_str(FALSE));
310 }
311 pr();
312 }
313
314
print_wiggle_map(map,base_like,scale)315 void print_wiggle_map(map,base_like,scale)
316 QTL_MAP *map;
317 real base_like, scale;
318 {
319 int i, j, stars;
320
321 if (map==NULL || scale<=0.0) send(CRASH);
322 i= map->num_intervals-1;
323
324 if (raw.data_type==BACKCROSS)
325 sf(ps,B1_WIGGLE_MAP,dist_str(map->qtl_pos[i],TRUE),
326 rsn(6.3,map->qtl_weight[i]),
327 map->var_explained*100.0,rsd(7.3,map->log_like));
328 else /* data_type==INTERCROSS */
329 sf(ps,F2_WIGGLE_MAP,dist_str(map->qtl_pos[i],TRUE),
330 rsn(6.3,map->qtl_weight[i]),rsn(6.3,map->qtl_dominance[i]),
331 map->var_explained*100.0,rsd(7.3,map->log_like));
332 pr();
333
334 if (map->log_like>base_like) {
335 stars=iminf((int)((map->log_like-base_like+scale+0.001)/scale),
336 (raw.data_type==BACKCROSS ? B1_SPACES_LEFT:F2_SPACES_LEFT)+1);
337 for (j=0; j<stars; j++) print("*");
338 }
339 nl();
340 }
341
342
print_wiggle_genetics(genetics)343 void print_wiggle_genetics(genetics)
344 GENETICS *genetics;
345 {
346 if (constrained(genetics)) {
347 print("Scanned QTL genetics are constrained to be: ");
348 print(genetics_str(genetics,TRUE)); nl();
349 } else print("Scanned QTL genetics are free.\n");
350 }
351
352
353 /* left_seq_str() in this file STILL needs CONT_VAR hooks! */
354
print_wiggle_left_seq(map)355 void print_wiggle_left_seq(map)
356 QTL_MAP *map;
357 {
358 char *left_seq;
359
360 if (nullstr(left_seq=left_seq_str(map))) print("No fixed-QTLs.\n");
361 else { print("Fixed-QTLs: "); print(left_seq); nl(); }
362 }
363
364
365 #define WIG_LIST_TITLE " NUM TRAIT%s SEQUENCE\n"
366 #define WIG_LIST_NAMES "%3d.%s %-10s %s\n"
367 #define WIG_LIST_NUMS "%3d.%s %3d %s\n"
368
print_saved_wiggles()369 void print_saved_wiggles()
370 {
371 int i;
372 WIGGLE_OPERATION *op;
373
374 if (first_wiggle==num_wiggles) return;
375 sf(ps,WIG_LIST_TITLE,(print_names ? " ":"")); pr();
376 for (i=first_wiggle; i<num_wiggles; i++) {
377 if ((op=wiggles[i])==NULL) send(CRASH);
378 if (op->data == NULL)
379 sf(ps," <deleted>\n");
380 else if (print_names)
381 sf(ps,WIG_LIST_NAMES,i+1,(op->num_orders==1 ? "1":"x"),
382 raw.trait_name[op->trait],op->seq_string);
383 else sf(ps,WIG_LIST_NUMS,i+1,(op->num_orders==1 ? "1":"x"),
384 op->trait+1,op->seq_string);
385 pr();
386 }
387 }
388
389
390 #define ORD_LIST_TITLE " NUM GENETICS FIXED-QTLS\n"
391 #define ORD_LIST_FORM "%3d.%-3d %-9s %s\n"
392
print_saved_wiggle(wiggle)393 void print_saved_wiggle(wiggle)
394 int wiggle;
395 {
396 int i, k;
397 WIGGLE_OPERATION *op;
398 char *left;
399
400 if (wiggle<0) return;
401 if (wiggle>=num_wiggles || (op=wiggles[wiggle])==NULL) send(CRASH);
402
403 print(ORD_LIST_TITLE);
404 k=op->num_intervals-1;
405 for (i=0; i<op->num_orders; i++) {
406 sf(ps,ORD_LIST_FORM,wiggle+1,i+1,
407 genetics_str(&op->data[i][0]->map->constraint[k],FALSE),
408 left_seq_str(op->data[i][0]->map));
409 pr();
410 }
411 }
412
413
print_saved_wiggle_order(wiggle,order,base_like,scale)414 void print_saved_wiggle_order(wiggle,order,base_like,scale)
415 int wiggle, order;
416 real base_like, scale;
417 {
418 int i, j, k, stars;
419 WIGGLE_OPERATION *op;
420 WIGGLE_INTERVAL **data;
421 WIGGLE_POINT *point;
422
423 if (wiggle<0) return;
424 if (wiggle>=num_wiggles || (op=wiggles[wiggle])==NULL) send(CRASH);
425 if (order<0 || order>=op->num_orders || base_like<0.0 || scale<=0.0)
426 send(CRASH);
427 data=op->data[order];
428
429 print_wiggle_title();
430 for (i=0; i<op->num_wiggled_intervals; i++) {
431 if (i>0 && !data[i]->contig) { print_wiggle_interval(NULL); nl(); }
432 print_wiggle_interval(data[i]->map);
433 for (j=0; j<data[i]->num_points; j++) {
434 point=data[i]->point[j];
435 if (raw.data_type==BACKCROSS)
436 sf(ps,B1_WIGGLE_MAP,dist_str(point->qtl_pos,TRUE),
437 rsn(6.3,point->qtl_weight),point->var_explained*100.0,
438 rsd(7.3,point->lod_score));
439 else /* data_type==INTERCROSS */
440 sf(ps,F2_WIGGLE_MAP,dist_str(point->qtl_pos,TRUE),
441 rsn(6.3,point->qtl_weight),rsn(6.3,point->qtl_dominance),
442 point->var_explained*100.0,rsd(7.3,point->lod_score));
443 pr();
444 if (point->lod_score>base_like) {
445 stars=
446 iminf((int)((point->lod_score-base_like+scale+0.001)/scale),
447 (raw.data_type==BACKCROSS ?B1_SPACES_LEFT:F2_SPACES_LEFT)+1);
448 for (k=0; k<stars; k++) print("*");
449 }
450 nl();
451 }
452 }
453 print_wiggle_interval(NULL);
454 }
455
456
457 #define PEAK_TITLE "QTL-Map for peak %d:\n"
458 #define PEAK_L_CONF "Confidence Interval: Left Boundary= "
459 #define PEAK_R_CONF " Right Boundary= "
460 #define PEAK_BOUND "%s + %s\n"
461 #define PEAK_NOBOUND "%s (off end)\n"
462
print_peak(peak,num)463 void print_peak(peak,num)
464 WIGGLE_PEAK *peak;
465 int num;
466 {
467 char *i;
468
469 print_map_divider();
470 sf(ps,PEAK_TITLE,num+1); pr();
471 print(PEAK_L_CONF);
472 if (peak->backward_pos==OFF_END) {
473 i=interval_str(peak->backward_left,peak->backward_right,FALSE);
474 sf(ps,PEAK_NOBOUND,i); pr();
475 } else {
476 i=interval_str(peak->backward_left,peak->backward_right,FALSE);
477 sf(ps,PEAK_BOUND,i,dist_str(peak->backward_pos,FALSE)); pr();
478 }
479 print(PEAK_R_CONF);
480 if (peak->forward_pos==OFF_END) {
481 i=interval_str(peak->forward_left,peak->forward_right,FALSE);
482 sf(ps,PEAK_NOBOUND,i); pr();
483 } else {
484 i=interval_str(peak->forward_left,peak->forward_right,FALSE);
485 sf(ps,PEAK_BOUND,i,dist_str(peak->forward_pos,FALSE)); pr();
486 }
487 nl();
488 if (peak->map==NULL) send(CRASH); /* KLUDGE for now */
489 print_qtl_map(peak->map,NULL);
490 }
491
492
493
494
495 /* SAVED TEST WIGGLES OUTPUT FORMAT... (F2 DATA ONLY)
496 -----------------------------------------------------------------------------|
497 GENETICS: FREE | DOMINANT | RECESSIVE | ADDITIVE |
498 -----------------------------------------------------------------------------|
499 LOG | LIKE | LIKE | LIKE |
500 POS WEIGHT DOM %VAR LIKE | %VAR DIFF | %VAR DIFF | %VAR DIFF |
501 -----------------------------------------------------------------------------|
502 interval= xxxxx length= 12.45 cm |
503 -----------------------------------------------------------------------------|
504 123.4 -2.345 -2.345 12.4 123.56 | 12.4 -12.67 | 12.4 -12.56 | 12.4 -12.56 |
505 -----------------------------------------------------------------------------|
506 */
507
508 #define TEST_WIGGLE_DIVIDER \
509 "-----------------------------------------------------------------------------|"
510 #define TEST_WIGGLE_TITLE \
511 "GENETICS: FREE | DOMINANT | RECESSIVE | ADDITIVE |"
512 #define TEST_WIGGLE_HEADER1 \
513 " LOG | LIKE | LIKE | LIKE |"
514 #define TEST_WIGGLE_HEADER2 \
515 "POS WEIGHT DOM %VAR LIKE | %VAR DIFF | %VAR DIFF | %VAR DIFF |"
516
517 #define TEST_WIGGLE_INTERVAL1 "interval= %s length= %s %s"
518 #define TEST_WIGGLE_MAP1 "%5s %s %s %4.1lf %s | "
519 #define TEST_WIGGLE_MAP2 "%4.1lf %s | %4.1lf %s | %4.1lf %s |%s"
520
521 /* We had to split up the TEST_WIGGLE_MAP line because of a (compiler?) bug
522 on the HP800. Leave it as is. */
523
524
print_test_wiggle_map(point,threshold)525 void print_test_wiggle_map(point,threshold)
526 WIGGLE_POINT **point;
527 real threshold;
528 {
529 real lod; lod= point[FREE]->lod_score;
530
531 sf(ps,TEST_WIGGLE_MAP1,
532 dist_str(point[FREE]->qtl_pos,TRUE),
533 rsn(6.3,point[FREE]->qtl_weight),
534 rsn(6.3,point[FREE]->qtl_dominance),
535 point[FREE]->var_explained*100.0,
536 rsd(6.2,point[FREE]->lod_score)); pr();
537
538 sf(ps,TEST_WIGGLE_MAP2,
539 point[DOMINANT]->var_explained*100.0,
540 rsd(6.2,point[DOMINANT]->lod_score-lod),
541 point[RECESSIVE]->var_explained*100.0,
542 rsd(6.2,point[RECESSIVE]->lod_score-lod),
543 point[ADDITIVE]->var_explained*100.0,
544 rsd(6.2,point[ADDITIVE]->lod_score-lod),
545 (lod+0.00499>=threshold ? "*\n":"\n")); pr();
546 }
547
548
print_test_wiggle_interval(map)549 void print_test_wiggle_interval(map)
550 QTL_MAP *map;
551 {
552 int i, last;
553 print(TEST_WIGGLE_DIVIDER); nl();
554 if (map==NULL) return;
555
556 last= map->num_intervals-1;
557 sf(ps,TEST_WIGGLE_INTERVAL1,
558 interval_str(map->left[last],map->right[last],FALSE),
559 dist_str(map->interval_len[last],FALSE),units_str(FALSE)); pr();
560 for (i=len(ps);i<77; i++) print(" ");
561 print("|\n"); print(TEST_WIGGLE_DIVIDER); nl();
562 }
563
564
print_test_wiggle_title()565 void print_test_wiggle_title()
566 { print(TEST_WIGGLE_DIVIDER); nl(); print(TEST_WIGGLE_TITLE); nl();
567 print(TEST_WIGGLE_DIVIDER); nl(); print(TEST_WIGGLE_HEADER1); nl();
568 print(TEST_WIGGLE_HEADER2); nl(); }
569
570
print_test_wiggle_order(wiggle,order,threshold)571 void print_test_wiggle_order(wiggle,order,threshold)
572 int wiggle, order;
573 real threshold;
574 {
575 int i, j, k, test;
576 WIGGLE_OPERATION *op;
577 WIGGLE_INTERVAL **data[NUM_MODELS];
578 WIGGLE_POINT *point[NUM_MODELS];
579
580 if (wiggle<0) return;
581 if (wiggle>=num_wiggles || (op=wiggles[wiggle])==NULL) send(CRASH);
582 if (order<0 || order>=op->num_orders || threshold<0.0)
583 send(CRASH);
584 for (test=0; test<NUM_MODELS; test++) data[test]=op->data[order+test];
585
586 print_test_wiggle_title();
587 for (i=0; i<op->num_wiggled_intervals; i++) {
588 if (i>0 && !data[FREE][i]->contig)
589 { print_test_wiggle_interval(NULL); nl(); }
590 print_test_wiggle_interval(data[FREE][i]->map);
591
592 for (j=0; j<data[FREE][i]->num_points; j++) {
593 for (test=0; test<NUM_MODELS; test++)
594 point[test]=data[test][i]->point[j];
595 print_test_wiggle_map(point,threshold);
596 }
597 }
598 print_test_wiggle_interval(NULL);
599 }
600
601
602
603 /*****************************************************************************/
604
print_trait(for_num_maps)605 void print_trait(for_num_maps)
606 int for_num_maps;
607 { sf(ps,"QTL map%s for trait %s:\n",maybe_s(for_num_maps),trait_str()); pr(); }
608
609
print_seq()610 void print_seq()
611 { print("Sequence: "); print(ints_string); nl(); }
612
613
print_old_seq(str)614 void print_old_seq(str)
615 char *str;
616 { print("Sequence: "); print(str); nl(); }
617
618
trait_str()619 char *trait_str()
620 {
621 char *str= get_temp_string();
622 sf(str,"%d (%s)",trait+1, raw.trait_name[trait]);
623 return(str);
624 }
625
626
interval_str(left,right,fill)627 char *interval_str(left,right,fill)
628 int left, right;
629 bool fill;
630 {
631 char *str= get_temp_string();
632 int i;
633 if (print_mapm_loci) {
634 if (!print_names) sf(str,"%d-%d",raw.original_locus[left],
635 raw.original_locus[right]);
636 else sf(str,"%s-%s",raw.locus_name[left],
637 raw.locus_name[right]);
638 if (fill) pad_to_len(str,print_names ? 21:9);
639 return(str);
640 }
641 else {
642 if (!print_names) sf(str,"%d-%d",left+1,right+1);
643 else sf(str,"%s-%s",raw.locus_name[left],raw.locus_name[right]);
644 if (fill) pad_to_len(str,print_names ? 21:9);
645 return(str);
646 }
647 }
648
dist_str(rec_frac,fill)649 char *dist_str(rec_frac,fill) /* Always takes 5 spaces */
650 real rec_frac;
651 bool fill; /* unused? */
652 {
653 char *str; str=get_temp_string();
654
655 if (units == CENTIMORGANS) {
656 if (fill) sf(str,"%-5.1lf",haldane_cm(rec_frac));
657 else sf(str,"%-3.1lf",haldane_cm(rec_frac));
658 str[5]='\0'; /* just in case */
659 } else sf(str,"%-5.3lf",rec_frac);
660
661 return(str);
662 }
663
664
units_str(fill)665 char *units_str(fill) /* Takes 5 spaces, if fill, 2 otherwise */
666 bool fill;
667 {
668 char *str; str=get_temp_string();
669
670 if (units == CENTIMORGANS) strcpy(str,"cM ");
671 else strcpy(str,"RF ");
672 if (!fill) str[2]='\0';
673
674 return(str);
675 }
676
677
genetics_str(genetics,verbose)678 char *genetics_str(genetics,verbose)
679 GENETICS *genetics;
680 bool verbose;
681 {
682 char *str= get_temp_string();
683
684 if (!data_loaded()) strcpy(str,"");
685
686 else if (raw.data_type==BACKCROSS) {
687 if (genetics->backx_weight==DONT_FIX) strcpy(str,"free");
688 else if (verbose) sf(str,"W=%-5.3lf",genetics->backx_weight);
689 else strcpy(str,"fixed");
690
691 } else if (raw.data_type==INTERCROSS)
692 switch (genetics->interx_type) {
693 case FREE: strcpy(str,"free"); break;
694 case DOMINANT: strcpy(str,"dominant"); break;
695 case ADDITIVE: strcpy(str,"additive"); break;
696 case RECESSIVE: strcpy(str,"recessive"); break;
697 case F3DOMINANT: strcpy(str,"f3dom"); break;
698 case F3RECESSIVE: strcpy(str,"f3rec"); break;
699 case TEST_MODELS: strcpy(str,"try"); break;
700 case CONSTRAINED:
701 if (!verbose) strcpy(str,"constrain");
702 else sf(str,"constraints: A=%-4.2lf B=%-4.2lf C=-4.2lf",
703 genetics->a,genetics->b,genetics->c);
704 break;
705 case FIXED:
706 if (!verbose) strcpy(str,"fixed");
707 else sf(str,"fixed: A=%-5.3lf D=%-5.3lf",genetics->a,
708 genetics->b);
709 break;
710 default: send(CRASH);
711 }
712
713 else send(CRASH);
714 return(str);
715 }
716
717
left_seq_str(map)718 char *left_seq_str(map)
719 QTL_MAP *map;
720 {
721 char *str=get_temp_string(), *qtl=get_temp_string();
722 char *interval, *genetics, pos[11];
723 int j;
724
725 for (j=0; j<map->num_intervals-1; j++) { /* all but rightmost interval */
726 interval= interval_str(map->left[j],map->right[j],FALSE);
727 genetics= genetics_str(&map->constraint[j],FALSE);
728 if (streq(genetics,"free")) genetics[0]='\0';
729 else strins(genetics," :");
730 if (map->fix_pos[j]==DONT_FIX) pos[0]='\0';
731 else sprintf(pos,"+%s",dist_str(map->fix_pos[j],FALSE));
732 sprintf(qtl,"%s[%s%s%s]",(str[0]!='\0' ?" ":""),interval,pos,genetics);
733 maxstrcat(str,qtl,TEMP_STRING_LEN);
734 }
735 return(str);
736 }
737
738
739
740 /********************* TINY QTL-MAP FORMAT: ****************************
741
742 INTERVALS %VAR LOG-LIKE
743 #123 13-56 23-56 33-56 43-56 53-56 63-56 73-56 83-56 93-56 12.4 123.56
744 POS: 00.0
745 W: 0.000
746 D: 0.000
747
748 #123 123-567 223-567 323-567 423-567 523-567 623-567 723-567 12.4 123.56
749 #123 1234567890-1234567890 1234567890-1234567890 12.4 1234.67 */
750
751 #define TINY_TITLE1 " INTERVAL%s:"
752 #define TINY_TITLE2 " %%VAR %s\n"
753 #define TINY_LINE1 "#%-3d "
754 #define TINY_LINE2 " %4.1lf %-7.2lf\n"
755
print_tiny_map(map,num,offset)756 void print_tiny_map(map,num,offset) /* unused, and broken */
757 QTL_MAP *map;
758 int num; /* number for map (starts at 0) or num_intervals if map is null */
759 real offset;
760 {
761 int i, spaces_per_int, ints_per_line, columns, right, n_ints;
762 char *interval;
763
764 if (map==NULL || map->num_continuous_vars>0) send(CRASH);
765
766 if (!print_names) {
767 spaces_per_int=raw.n_loci<=99 ? 6:8;
768 ints_per_line= raw.n_loci<=99 ? 9:7;
769 } else { spaces_per_int=22; ints_per_line=2; }
770
771 columns= (spaces_per_int * map->num_intervals);
772 if (columns>60) columns=60;
773 right=imaxf(5+columns,15);
774
775 if (num<0) { /* header only */
776 sf(ps,TINY_TITLE1,n_ints>1 ? "S":" "); pr(); to_column(right);
777 sf(ps,TINY_TITLE2,(offset>0.0 ? "LIKE-DIFF":"LOG-LIKE")); pr();
778
779 } else { /* print the map */
780 sf(ps,TINY_LINE1,num+1); pr();
781 for (i=0; i<map->num_intervals; i++) {
782 interval= interval_str(map->left[i],map->right[i],FALSE);
783 print(interval); print(" ");
784 if ((i+1)%ints_per_line==0 && i!=map->num_intervals-1)
785 print("\n ");
786 }
787 to_column(right);
788 sf(ps,TINY_LINE2,map->var_explained*100.0,map->log_like-offset); pr();
789 }
790 }
791
792
793
794 /***** Show-compares printing functions *****/
795
796 #define COMP_LIST_TITLE " NUM TRAIT%s SEQUENCE\n"
797 #define COMP_LIST_NAMES "%3d.%s %-10s %s\n"
798 #define COMP_LIST_NUMS "%3d.%s %3d %s\n"
799
print_saved_compares()800 void print_saved_compares()
801 {
802 int i;
803 COMPARE_OPERATION *op;
804
805 if(first_compare == num_compares) return;
806 sf(ps,COMP_LIST_TITLE,(print_names ? " ":"")); pr();
807 for(i=first_compare; i<num_compares; i++) {
808 if ((op=compares[i]) == NULL) send(CRASH);
809 if (op->data==NULL)
810 sf(ps," This compare (number%d) has been deleted!\n",i+1);
811 else if(print_names)
812 sf(ps,COMP_LIST_NAMES,i+1,(op->num_contigs==1 ? "1" : "x"),
813 raw.trait_name[op->trait],op->seq_string);
814 else
815 sf(ps,COMP_LIST_NUMS,i+1,(op->num_contigs==1 ? "1" : "x"),
816 op->trait+1,op->seq_string);
817 pr();
818 }
819 }
820
821
822 #define GOOD_MAPS_TITLE "QTL-maps above LOD %sfalloff:\n"
823 #define BAD_MAPS_TITLE "QTL-maps below LOD %sfalloff:\n"
824 #define COMP_MAP_TOP "QTL-map #%d: LOD score difference= %-5.2lf\n"
825
print_best_saved_maps(compare,contig,threshold,falloff)826 void print_best_saved_maps(compare,contig,threshold,falloff)
827 int compare, contig;
828 real threshold, falloff;
829 {
830 COMPARE_OPERATION *op;
831 int start, i, *like_index, num_in_list, j;
832 bool bad_maps_yet, good_maps_printed;
833 real best, *like_list;
834
835 like_index= NULL; like_list= NULL;
836 run {
837 best= 0.0;
838 array(like_list,raw.n_loci,real);
839 array(like_index,raw.n_loci,int);
840
841 if ((op=compares[compare])==NULL) send(CRASH);
842
843 for (i=0, start=0; i<contig; i++)
844 { while (op->data[start]->contig) start++; }
845
846 like_list[0]=op->data[start]->map->log_like;
847 like_index[0]=start;
848 num_in_list=1;
849 start++;
850
851 /* while (op->data[start]->contig) { */
852
853 for(j = start; j<op->num_orders; j++) {
854 i=num_in_list;
855 while (op->data[start]->map->log_like > like_list[i-1]) {
856 like_list[i]=like_list[i-1];
857 like_index[i]=like_index[i-1];
858 i--;
859 if (i==0) break;
860 }
861 like_list[i]=op->data[start]->map->log_like;
862 like_index[i]=start;
863 num_in_list++; start++;
864 }
865 best=like_list[0];
866
867 bad_maps_yet=good_maps_printed=FALSE;
868 if (best>=threshold) {
869 sf(ps,GOOD_MAPS_TITLE,(falloff<0.0 ? "threshold minus ":"")); pr();
870 print_map_divider();
871 good_maps_printed=TRUE;
872 }
873
874 for (i=0; i<num_in_list; i++) {
875 map=op->data[like_index[i]]->map;
876 if (good_maps_printed && ((falloff>=0.0 && map->log_like>falloff)||
877 (falloff<0.0 && map->log_like > best+falloff))) {
878 sf(ps,COMP_MAP_TOP,i+1,map->log_like-best); pr(); nl();
879 print_qtl_map(map,NULL);
880 print_map_divider();
881 } else {
882 if (!bad_maps_yet) {
883 if (good_maps_printed) nl();
884 sf(ps,BAD_MAPS_TITLE,
885 (falloff<0.0 ? "threshold minus ":"")); pr();
886 print_tiny_map(map,-1,0.0); /* title */
887 bad_maps_yet=TRUE;
888 }
889 print_tiny_map(map,i,0.0);
890 }
891 }
892
893 } on_exit {
894 unarray(like_index,int);
895 unarray(like_list,real);
896 relay_messages;
897 }
898 }
899
900
901 #define SAVED_BEST_LIKE \
902 "LOD score maximum for these maps: %-4.2lf\n"
903
print_saved_maps(compare,contig)904 void print_saved_maps(compare,contig)
905 int compare, contig;
906 {
907 COMPARE_OPERATION *op;
908 int i, start;
909 bool bad_maps_yet, good_maps_printed;
910 real best, like;
911
912 if ((op=compares[compare])==NULL) send(CRASH);
913 best= 0.0;
914
915 for (i=0, start=0; i<contig; i++)
916 { while (op->data[start]->contig) start++; }
917 for (i=start; i<op->num_orders; i++) {
918 if (i>start && !op->data[i]->contig) break;
919 if ((like=op->data[i]->map->log_like)>best) best=like;
920 }
921
922 sf(ps,SAVED_BEST_LIKE,best); pr(); nl();
923 print_tiny_map(op->data[0]->map,-1,1.0); /* title */
924
925 for (i=start; i<op->num_orders; i++) {
926 if (i>start && !op->data[i]->contig) break;
927 print_tiny_map(op->data[i]->map,i,best);
928 }
929 }
930
get_fixed_qtl_weights(map)931 void get_fixed_qtl_weights(map)
932 QTL_MAP *map;
933 {
934 int i;
935
936 if (raw.data_type==BACKCROSS) print("Enter weight for QTL(s):\n");
937 else print("Enter weight and dominance for QTL(s):\n");
938 for (i=0; i<map->num_intervals; i++) {
939 print(interval_str(map->left[i],map->right[i],FALSE));
940 getln(": ");
941 if (!rtoken(&ln,rREQUIRED,&map->qtl_weight[i])) error("bad weight!");
942 if (raw.data_type==INTERCROSS &&
943 !rtoken(&ln,rREQUIRED,&map->qtl_dominance[i]))error("bad weight!");
944 }
945 print("Enter weight for cont-vars:\n");
946 for (i=0; i<map->num_continuous_vars; i++) {
947 if (map->cont_var[i]==EPISTASIS_TERM) strcpy(ps,"epistasis: ");
948 else sf(ps,"%d (%s): ",map->cont_var[i]+1,raw.trait_name[map->cont_var[i]]);
949 getln(ps);
950 if (!rtoken(&ln,rREQUIRED,&map->cont_var_weight[i]))
951 error("bad weight!");
952 }
953 /* getln("Enter mu and sigma-squared: ");
954 if (!rtoken(&ln,rREQUIRED,&map->mu) ||
955 !rtoken(&ln,rREQUIRED,&map->sigma_sq)) error("bad weight!"); */
956 }
957
958
959
960
961