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_EQN
17 #define INC_STATS
18 #define INC_CALLQCTM
19 #define INC_QTOPLEVEL
20 #define INC_QLOWLEVEL
21 #include "qtl.h"
22
23 /***** GLOBAL *****/
24 int print_long_maps;
25 char default_intercross_chars[10], default_backcross_chars[10];
26 char geno_chars[10];
27 bool auto_save;
28 bool fix_weight_kludge;
29 real min_trait_val, max_trait_val;
30
31 /***** INTERNAL *****/
32 char *temp, errmsg[MAXLINE+1];
33 int get_pair_entry();
34 void load_qtl_files(), save_qtl_files();
35 bool save_on_exit();
36 void print_ps_wiggle_order(), print_ps_multi_wiggle();
37
38 /***** QTL Commands in this file *****/
39 command set_intervals();
40 command show_seq_history();
41 command set_trait();
42 command set_genetics();
43 command load_data();
44 command qtl_map();
45 command singles(); /* FROB */
46 command qtl_like();
47 command compare();
48 command pheno_histogram();
49 command wiggle();
50 command test_wiggle();
51 command make_trait();
52 command predict();
53 command forget_trait();
54 command forget_compare();
55 command forget_all_compares();
56 command forget_all_scans();
57 command forget_scan();
58 command set_units();
59 command set_print_iter();
60 command set_print_long();
61 command set_dbg_qctm();
62 command set_max_ints();
63 command set_brute();
64 command set_print_bf();
65 command set_print_nam();
66 command set_print_scans();
67 command set_print_rm();
68 command set_tolerance();
69 command set_bag_qctm();
70 command set_autosave();
71 command set_more_mode();
72 command translate();
73 command let();
74 command names();
75 command forget();
76 command sequence_editor();
77 command show_peaks();
78 command show_wiggle();
79 command draw_wiggle();
80 command show_trait();
81 command list_traits();
82 command list_wiggles();
83 command save_status();
84 command show_trait();
85 command show_test_wiggle();
86 command show_best();
87 command show_compare();
88 command list_compares();
89 command new_show_map();
90 command tester();
91 command dump_traits();
92 command dump_genome();
93 command dump_scan();
94 command tweak_weight();
95 command set_min_trait();
96 command set_max_trait();
97
98 /* QTL Command Topics */
99 #define RAWDATA 1
100 #define SELECT 2
101 #define MAPPING 3
102 #define DATAVIEW 4
103 #define OPTION 5
104 #define SYSTEM 6
105 #define ALIAS 7
106 #define WIZARD 8
107 #define HELPONLY 9
108
109 DATA *data;
110
111 /***** Initialize all user variables and the command tables *****/
cmd_init()112 void cmd_init() {
113
114 array(temp,MAXLINE+1,char);
115
116 /********** NO LONGER USE THIS ***********************************************
117 mktopic(SELECT, "Interval (sequence) and Trait Selection Command", CMD);
118 mktopic(MAPPING, "QTL-Map Making Function", CMD);
119 mktopic(DATAVIEW,"QTL-Map Data Viewing Function", CMD);
120 mktopic(OPTION, "MAPMAKER/QTL Option", OPT);
121 mktopic(SYSTEM, "System Command", CMD);
122 mktopic(RAWDATA, "Raw Data Handling Command", CMD);
123 mktopic(ALIAS, "Command Aliase", WIZ);
124 mktopic(WIZARD, "Wizard Command", WIZ);
125 mktopic(HELPONLY,"Other Help Topics", HLP);
126 *****************************************************************************/
127
128 /* COMMANDS: command-name abbrev c-function topic code*/
129 /* 1234567890123456789012345 123 1234567890123456 12345678 123*/
130 /* --------------------------- ----- ---------------- -------- ---*/
131 /* mkcommand("1234567890123456789012345","123",1234567890123456,12345678,1);*/
132
133 /* system commands in shell.c... */
134 mkcommand("help", "?", help, CMD);
135 mkcommand("quit", "q", quit, CMD);
136 mkcommand("run", "r", run_from_file, CMD);
137 mkcommand("photo", "", do_photo, CMD);
138 mkcommand("change directory", "cd", cd_command, CMD);
139 mkcommand("system", "!", system_command, CMD);
140 mkcommand("previous commands", "", show_cmd_history,CMD);
141 mkcommand("review output", "", review_output, CMD);
142 mkcommand("time", "", show_time, CMD);
143 mkcommand("wizard mode", "", set_wizard, CMD);
144 mkcommand("comment", "", comment, CMD);
145 mkcommand("about mapmaker/qtl", "", about, CMD);
146
147 /* QTL specific commands in this file... */
148 mkcommand("scan", "", wiggle, CMD);
149 mkcommand("map", "m", qtl_map, CMD);
150 mkcommand("like", "", qtl_like, WIZ);
151 mkcommand("compare", "", compare, CMD);
152
153 mkcommand("list scans", "", list_wiggles, CMD);
154 mkcommand("show scan", "", show_wiggle, CMD);
155 mkcommand("draw scan", "", draw_wiggle, CMD);
156 mkcommand("show peaks", "", show_peaks, CMD);
157 mkcommand("show trys", "", show_test_wiggle,CMD);
158 mkcommand("forget scan", "", forget_scan, CMD);
159 mkcommand("forget all scans", "", forget_all_scans,CMD);
160
161 mkcommand("list compares", "", list_compares, CMD);
162 mkcommand("show best maps", "", show_best, CMD);
163 mkcommand("show compare", "", show_compare, CMD);
164
165 /***** ALIASES - ELIMINATE
166 mkcommand("wiggle", "", wiggle, CMD);
167 mkcommand("list wiggles", "", list_wiggles, CMD);
168 mkcommand("show wiggle", "", show_wiggle, CMD);
169 mkcommand("draw wiggle", "", draw_wiggle, CMD);
170 *****/
171
172 mkcommand("sequence", "s", set_intervals, CMD);
173 mkcommand("history", "h", show_seq_history,CMD);
174 mkcommand("trait", "t", set_trait, CMD);
175 mkcommand("let", "", let, CMD);
176 mkcommand("names", "", names, CMD);
177 mkcommand("forget named sequence", "", forget, CMD);
178 mkcommand("edit sequence", "", sequence_editor, CMD);
179
180
181 /* mkcommand("prepare data", "", prep_data, CMD); */
182 mkcommand("load data", "ld", load_data, CMD);
183 mkcommand("save data", "", save_status, CMD);
184 mkcommand("auto save data", "", set_autosave, OPT);
185 mkcommand("list loci", "", translate, CMD);
186 mkcommand("list traits", "", list_traits, CMD);
187 mkcommand("show linkage map", "", new_show_map, CMD);
188 mkcommand("show trait", "", show_trait, CMD);
189 mkcommand("make trait", "", make_trait, CMD);
190 mkcommand("forget trait", "", forget_trait, CMD);
191
192 mkcommand("forget compare", "", forget_compare, CMD);
193 mkcommand("forget all compares", "",forget_all_compares,CMD);
194
195 mkcommand("units", "", set_units, PAR);
196 mkcommand("print names", "", set_print_nam, OPT);
197 mkcommand("print scans", "", set_print_scans, WIZ);
198 mkcommand("more mode", "", set_more_mode, OPT);
199
200 mkcommand("print rec mat", "", set_print_rm, WIZ);
201 mkcommand("brute force", "", set_brute, WIZ);
202 mkcommand("bag qctm", "", set_bag_qctm, WIZ);
203 mkcommand("debug qctm", "", set_dbg_qctm, WIZ);
204 mkcommand("print brute force", "", set_print_bf, WIZ);
205 mkcommand("print bf", "", set_print_bf, WIZ);
206 mkcommand("print long maps", "plm",set_print_long, WIZ);
207 mkcommand("print iterations", "pi", set_print_iter, WIZ);
208 mkcommand("tolerance", "", set_tolerance, WIZ);
209
210 mkcommand("test", "", tester, WIZ);
211 mkcommand("predict trait", "", predict, WIZ);
212 mkcommand("dump traits", "", dump_traits, WIZ);
213 mkcommand("dump genome", "", dump_genome, WIZ);
214 mkcommand("dump scan", "", dump_scan, WIZ);
215 mkcommand("tweak weight", "", tweak_weight, WIZ);
216 mkcommand("min trait", "", set_min_trait, WIZ);
217 mkcommand("max trait", "", set_max_trait, WIZ);
218
219 /* FROB */
220 mkcommand("singles", "", singles, WIZ);
221
222 /* QTL User state variables */ /* source file defined in */
223 units= CENTIMORGANS; /* qtop.c */
224 print_names= FALSE; /* qtop.c */
225 print_mapm_loci = TRUE; /* qtop.c */
226 print_scans = TRUE; /* qtop.c */
227 print_maps= FALSE; /* qtop.c */
228 print_long_maps= FALSE; /* qcmds.c - here */
229 trait= NOTRAIT; /* qcmds.c */
230 auto_save= TRUE; /* qcmds.c */
231 more_mode= TRUE; /* qcmds.c */
232
233 print_iter= FALSE; /* qctm.c */
234 print_rec_mat= FALSE; /* qctm.c */
235 debug_qctm= FALSE; /* qctm.c */
236 bag_qctm= FALSE; /* qctm.c */
237 pos_tolerance= 0.001; /* qctm.c */
238 like_tolerance= 0.001; /* qctm.c */
239 mat_tolerance= 0.00001; /* qctm.c */
240 brute_force= TRUE; /* qctm.c */
241 print_brute_force= FALSE; /* qctm.c */
242
243 /* ints, ints_string, and seq_history are set up in seq_init() */
244
245 /* Other global state variables */ /* source file */
246 max_intervals= -5; /* qctm.c *** Negative -> default val*/
247 max_continuous_vars= -5; /* ditto */
248 wizard_mode= FALSE; /* shell.c */
249 quit_save_hook= save_on_exit; /* shell.c - function to call */
250 fix_weight_kludge= FALSE;
251 min_trait_val= -VERY_BIG; max_trait_val=VERY_BIG;
252 }
253
254
help_init()255 void help_init()
256 {
257 /* #include "qtl.code" - outdated help mechanism */
258 }
259
260
261 /********** These are the QTL commands ***********/
262
set_print_iter()263 command set_print_iter() { maybe_set_bool(&print_iter); }
set_print_rm()264 command set_print_rm() { maybe_set_bool(&print_rec_mat); }
set_print_long()265 command set_print_long() { maybe_set_bool(&print_long_maps); }
set_print_nam()266 command set_print_nam() { maybe_set_bool(&print_names); }
set_print_scans()267 command set_print_scans(){ maybe_set_bool(&print_scans); }
set_brute()268 command set_brute() { maybe_set_bool(&brute_force); }
set_print_bf()269 command set_print_bf() { maybe_set_bool(&print_brute_force); }
set_dbg_qctm()270 command set_dbg_qctm() { maybe_set_bool(&debug_qctm); }
set_bag_qctm()271 command set_bag_qctm() { maybe_set_bool(&bag_qctm); }
set_tolerance()272 command set_tolerance() { maybe_set_real(&like_tolerance,0.0,VERY_BIG,8.6); }
set_autosave()273 command set_autosave() { maybe_set_bool(&auto_save); }
set_more_mode()274 command set_more_mode() { maybe_set_bool(&more_mode); }
275
set_min_trait()276 command set_min_trait() { maybe_set_real(&min_trait_val,0.0,VERY_BIG,8.6); }
set_max_trait()277 command set_max_trait() { maybe_set_real(&max_trait_val,0.0,VERY_BIG,8.6); }
278
set_max_ints()279 command set_max_ints() /* BUG? e can only do this before qctm_init! */
280 { maybe_set_int(&max_intervals,1,7); } /* limit is discussed in qctm.c */
281
set_units()282 command set_units()
283 {
284 if (nullstr(args)) {
285 if(units == RECFRACS)
286 print("the 'units' are currently set to recombination-fractions.\n");
287 else
288 print("the 'units' are currently set to haldane centimorgans.\n");
289
290 } else { /* have args */
291 crunch(args);
292 if (matches(args,"recombination fractions") ||
293 matches(args,"rec-fracs") ||
294 matches(args,"recombination-fractions") ||
295 matches(args,"rec fracs") ||
296 matches(args,"rf")) units = RECFRACS;
297 else if (matches(args,"centimorgans") ||
298 matches(args,"cm") ||
299 matches(args,"haldane centimorgans") ||
300 matches(args,"haldane cm")) units= CENTIMORGANS;
301 else usage_error(1);
302 if(units == RECFRACS)
303 print("the 'units' are now set to recombination-fractions.\n");
304 else
305 print("the 'units' are now set to haldane centimorgans.\n");
306 }
307 }
308
309
set_intervals()310 command set_intervals()
311 {
312 int errpos;
313
314 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
315 if (nullstr(args)) {
316 sf(ps,"The sequence is%s%s'\n",
317 (len(ints_string)<55 ? " '":":\n'"),ints_string); pr();
318
319 } else if (set_qtl_sequence(args,errmsg,&errpos)) {
320 sf(ps,"The sequence is now%s%s'",
321 (len(ints_string)<55 ? " '":":\n'"),ints_string);
322 maybe_ok(ps);
323
324 } else {
325 print("Error in sequence '"); print(args); print("'\n ");
326 space((errpos+20) % LINE); print("^\n"); print(errmsg); nl();
327 }
328 }
329
330
show_seq_history()331 command show_seq_history()
332 {
333 int i, num_to_print, first_to_print, printed_any;
334 char *cmd_str;
335 TABLE *seq_history;
336
337 get_one_arg(itoken,1000,&num_to_print);
338 if (num_to_print<1) usage_error(1);
339
340 printed_any=FALSE;
341 first_to_print= imaxf(context[active_context]->seq_history_num-num_to_print,0);
342 seq_history= context[active_context]->sequence_history;
343 for(Te=seq_history->list; Te!=NULL; Te=Te->next) {
344 i=Te->id.num; cmd_str=Te->string;
345 if (i>=first_to_print && i<context[active_context]->seq_history_num) {
346 if(!printed_any){print("Previous sequences:\n"); printed_any=TRUE;}
347 sf(ps,"%3d ",i+1); pr(); print(cmd_str); nl();
348 }
349 }
350 if (!printed_any) print("No sequences have yet been entered.\n");
351 }
352
353
set_trait()354 command set_trait()
355 {
356 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
357 run {
358 if (nullstr(args)) {
359 if (trait<0) print("The current trait is not set.\n");
360 else { print("The current trait is: "); print(trait_str()); nl(); }
361
362 } else { /* there was an argument */
363 set_trait_spec(args);
364 sf(ps,"The current trait is now: %s",trait_str());
365 maybe_ok(ps);
366 }
367 }
368 except_when(BADTRAIT) {
369 sf(ps,"Bad trait name or number specified.\n%s",BADTRAIT_errmsg);
370 error(ps);
371 }
372 }
373
374
translate()375 command translate()
376 {
377 int printed, i, j;
378
379 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
380 nomore_args(0);
381
382 printed= 0;
383 print(BIG_DIVIDER); print("LOCI:\n");
384 for (i=0; i<raw.n_loci; i++) {
385 if ((printed++)%5==0) nl();
386 sf(ps,"%3d %-10s ",raw.original_locus[i],raw.locus_name[i]); pr();
387 }
388 nl();
389
390 /* nl(); print("\nTRAITS:\n"); printed=0;
391 for (i=0; i<raw.n_traits; i++) {
392 if ((printed++)%5==0) nl();
393 if (nullstr(raw.trait_name[i])) sf(ps,"%3d <deleted> ",i+1);
394 else sf(ps,"%3d %-10s ",i+1,raw.trait_name[i]);
395 pr();
396 } nl(); */
397 print(BIG_DIVIDER);
398
399 /* Old code - what to do if there are args?
400 while (stoken(&args,sREQUIRED,temp)) {
401 if (valid_locus_str(temp,&i,&errmsg)) {
402 sf(ps," locus %d %-10s\n",i+1,raw.locus_name[i]); print(ps);
403 } else if (valid_trait_str(temp,&i,&errmsg)) {
404 sf(ps," trait %d %-10s\n",i+1,raw.trait_name[i]); print(ps);
405 } else {
406 sf(ps," Error: \"%s\" doesn't match, or is ambiguous.\n",
407 temp); print(ps);
408 }
409 } */
410 }
411
412
qtl_map()413 command qtl_map()
414 {
415 int perm;
416
417 qtl_ready(ANY_DATA,SEQ,TRAIT,QCTM);
418 nl();
419
420 for_all_orders(ints,map,perm) {
421 make_qtl_map(map);
422 print_map_divider();
423 print_trait(1); nl();
424 print_qtl_map(map,free_genetics);
425 }
426 print_map_divider();
427 }
428
429
singles()430 command singles() /* FROB */
431 {
432 int perm, last;
433 real threshold, scale;
434
435 qtl_ready(ANY_DATA,SEQ,TRAIT,QCTM);
436 get_arg(rtoken,2.0,&threshold);
437 get_arg(rtoken,0.25,&scale);
438 last= num_intervals-1;
439
440 nl();
441 print_trait(num_orders); nl();
442 nl();
443 print_short_title();
444
445 for_all_orders(ints,map,perm) {
446 map->fix_pos[last]= 0.0;
447 make_qtl_map(map);
448 print_short_qtl_map(map,threshold,scale);
449 }
450 nl();
451 }
452
453
qtl_like()454 command qtl_like()
455 {
456 int perm;
457 real like;
458
459 /* qtl_ready(ANY_DATA,SEQ,TRAIT,QCTM);
460 nl();
461
462 for_all_orders(ints,map,perm) {
463 like= qtl_map_like(map);
464 sf(ps,"like= %-7.2lf\n",like); pr();
465 } */
466 }
467
468
compare()469 command compare()
470 {
471 int perm, i, comp;
472
473 qtl_ready(ANY_DATA,SEQ,TRAIT,QCTM);
474 run {
475 comp=allocate_compare_struct(trait,ints,ints_string,num_intervals,
476 num_orders);
477
478 for_all_orders(ints,map,perm) {
479 if (perm==FIRST || perm>CONTIG)
480 { i=0; nl(); print_tiny_map(map,-1,0.0); } /* print tiny title */
481
482 make_qtl_map(map);
483 print_tiny_map(map,i++,0.0);
484 store_compare_map(comp,map,(perm<=CONTIG));
485 }
486 } on_exit {
487 if(msg == NOMEMORY || msg == INTERRUPT || msg == CRASH)
488 bash_compare_struct(num_compares-1);
489 }
490 }
491
492
list_compares()493 command list_compares()
494 {
495 if (first_compare==num_compares)
496 print("No compare results have been saved.\n");
497 else { print("\nSaved compare results:\n\n"); print_saved_compares(); }
498 }
499
500 #define SAVED_COMP_BEST "Sorted compare results number%s %d.%d"
501 #define SAVED_COMP_COMPS "-%d.%d"
502 #define SAVED_COMP_TRAIT " for trait %d (%s).\n"
503 #define SAVED_COMP_NUM "Saved compare number %d.%d:\n"
504 #define SAVED_COMP_THRESH "LOD Threshold: %-4.2lf Falloff: %-4.2lf\n"
505
show_best()506 command show_best()
507 {
508 char arg[TOKLEN+1];
509 int compare, contig, i;
510 real threshold, falloff;
511 COMPARE_OPERATION *op;
512
513 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
514 get_arg(stoken,"",arg);
515 get_arg(rtoken,2.0,&threshold);
516 get_arg(rtoken,-1.0,&falloff);
517 nomore_args(num_args);
518 if (threshold<0.0 || falloff==0.0 || falloff>=threshold)
519 usage_error(num_args);
520
521 get_compare_nums(arg,&compare,&contig);
522 if (compare<0) {
523 compare= num_compares-1;
524 contig= compares[compare]->num_contigs-1;
525 } else if (contig<0 && compares[compare]->num_contigs==1) contig=0;
526 op=compares[compare]; nl();
527
528 if (contig<0) {
529 sf(ps,SAVED_COMP_BEST,maybe_s(op->num_contigs),compare+1,1); pr();
530 if (op->num_contigs>1)
531 { sf(ps,SAVED_COMP_COMPS,compare+1,op->num_contigs); pr(); }
532 sf(ps,SAVED_COMP_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
533 print_old_seq(op->seq_string);
534 sf(ps,SAVED_COMP_THRESH,threshold,falloff); pr();
535 for (i=0; i<op->num_contigs; i++) {
536 nl(); sf(ps,SAVED_COMP_NUM,compare+1,i+1); pr();
537 print_best_saved_maps(compare,i,threshold,falloff);
538 }
539
540 } else {
541 sf(ps,SAVED_COMP_BEST,"",compare+1,contig+1); pr();
542 sf(ps,SAVED_COMP_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
543 print_old_seq(op->seq_string);
544 sf(ps,SAVED_COMP_THRESH,threshold,falloff); pr();
545 nl(); print_best_saved_maps(compare,contig,threshold,falloff);
546 }
547 }
548
549
550 #define SAVED_COMP_WHICH "Compare results number%s %d.%d"
551
show_compare()552 command show_compare()
553 {
554 char arg[TOKLEN+1];
555 int compare, contig, i;
556 COMPARE_OPERATION *op;
557
558 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
559 get_arg(stoken,"",arg);
560 nomore_args(num_args);
561
562 get_compare_nums(arg,&compare,&contig);
563 if (compare<0) {
564 compare= num_compares-1;
565 contig= compares[compare]->num_contigs-1;
566 } else if (contig<0 && compares[compare]->num_contigs==1) contig=0;
567 op=compares[compare]; nl();
568
569 if (contig<0) {
570 sf(ps,SAVED_COMP_WHICH,maybe_s(op->num_contigs),compare+1,1); pr();
571 if (op->num_contigs>1)
572 { sf(ps,SAVED_COMP_COMPS,compare+1,op->num_contigs); pr(); }
573 sf(ps,SAVED_COMP_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
574 print_old_seq(op->seq_string);
575 for (i=0; i<op->num_contigs; i++) {
576 nl(); sf(ps,SAVED_COMP_NUM,compare+1,i+1); pr();
577 print_saved_maps(compare,i);
578 }
579
580 } else {
581 sf(ps,SAVED_COMP_WHICH,"",compare+1,contig+1); pr();
582 sf(ps,SAVED_COMP_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
583 print_old_seq(op->seq_string);
584 nl(); print_saved_maps(compare,contig);
585 }
586 }
587
588
589 #define WIGGLE_STORED "Results have been stored as scan number %d.\n"
590 #define WIGGLES_STORED "Results have been stored as scan numbers %d.1-%d.%d.\n"
591 #define THRESHOLD_AND_SCALE "LOD threshold: %-4.2lf Scale: %-4.2lf per '*'\n"
592
wiggle()593 command wiggle()
594 {
595 int perm, wiggle;
596 real inc, threshold, scale;
597
598
599 qtl_ready(ANY_DATA,WIGSEQ,TRAIT,QCTM);
600 get_arg(rtoken,2.0,&inc);
601 get_arg(rtoken,2.0,&threshold);
602 get_arg(rtoken,0.25,&scale);
603 nomore_args(num_args);
604 if (inc<0.0 || threshold<0.0 || scale<0.0) usage_error(num_args);
605
606 if (print_scans) {
607 nl(); print_trait(2); print_seq();
608 sf(ps,THRESHOLD_AND_SCALE,threshold,scale); pr();
609 }
610 run {
611 wiggle=allocate_wiggle_struct(trait,ints,ints_string,num_intervals,
612 num_orders,num_ints_to_wiggle);
613
614 for_wiggle_orders(ints,map,inc,perm) {
615
616 if (print_scans) {
617 if (perm==FIRST || perm==TEST || perm==ORDER) {
618 if (perm!=FIRST) { print_wiggle_interval(NULL); }
619 nl(); print_wiggle_left_seq(map);
620 print_wiggle_genetics(&map->constraint[num_intervals-1]);
621 nl();
622 print_wiggle_title();
623 print_wiggle_interval(map);
624 store_wiggle_interval(wiggle,map,TRUE,FALSE,inc);
625 }
626 else if (perm==CONTIG || perm==SKIP) {
627 if (perm==SKIP) { print_wiggle_interval(NULL); nl(); }
628 print_wiggle_interval(map);
629 store_wiggle_interval(wiggle,map,FALSE,(perm==CONTIG),inc);
630 }
631 }
632 else {
633 if (perm==FIRST || perm==TEST || perm==ORDER)
634 store_wiggle_interval(wiggle,map,TRUE,FALSE,inc);
635 else if (perm==CONTIG || perm==SKIP)
636 store_wiggle_interval(wiggle,map,FALSE,(perm==CONTIG),inc);
637 }
638 make_qtl_map(map);
639 if (print_scans) print_wiggle_map(map,threshold,scale);
640 store_wiggle_point(wiggle,map);
641 }
642
643 if (print_scans)
644 print_wiggle_interval(NULL); nl(); /* line on the bottom */
645 if (wiggle>=0) {
646 if (num_orders==1) sf(ps,WIGGLE_STORED,wiggle+1);
647 else sf(ps,WIGGLES_STORED,wiggle+1,wiggle+1,num_orders);
648 pr();
649 }
650 } on_exit {
651 if(msg == NOMEMORY || msg == INTERRUPT || msg == CRASH)
652 bash_wiggle_struct(num_wiggles-1); /* the current one */
653 }
654 }
655
656
657
658 #define SAVED_WIGGLE_WIG "Scan results number%s %d.%d"
659 #define SAVED_WIGGLE_WIGS "-%d.%d"
660 #define SAVED_WIGGLE_TRAIT " for trait %d (%s).\n"
661
list_wiggles()662 command list_wiggles()
663 {
664 char arg[TOKLEN+1];
665 int wiggle, order, last;
666 WIGGLE_OPERATION *op;
667
668 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
669 get_arg(stoken,"",arg); nomore_args(num_args);
670
671 if (first_wiggle==num_wiggles) print("No scan results have been saved.\n");
672 if (nullstr(arg))
673 { print("\nSaved scan results:\n\n"); print_saved_wiggles(); return; }
674
675 get_wiggle_nums(arg,&wiggle,&order);
676 if (wiggle<0) { wiggle=num_wiggles-1; order=wiggles[wiggle]->num_orders-1;}
677 else if (order<0 && wiggles[wiggle]->num_orders==1) order=0;
678 op=wiggles[wiggle]; nl();
679
680 sf(ps,SAVED_WIGGLE_WIG,maybe_s(op->num_orders),wiggle+1,1); pr();
681 if (op->num_orders>1)
682 { sf(ps,SAVED_WIGGLE_WIGS,wiggle+1,op->num_orders); pr(); }
683 sf(ps,SAVED_WIGGLE_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
684 print_old_seq(op->seq_string); nl();
685 print_saved_wiggle(wiggle);
686 }
687
688
show_wiggle()689 command show_wiggle()
690 {
691 char arg[TOKLEN+1];
692 int wiggle, order, last;
693 real threshold, scale;
694 WIGGLE_OPERATION *op;
695
696 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
697 get_arg(stoken,"",arg);
698 get_arg(rtoken,2.0,&threshold);
699 get_arg(rtoken,0.25,&scale);
700 nomore_args(num_args);
701
702 if (threshold<0.0 || scale<=0.0) usage_error(2);
703 get_wiggle_nums(arg,&wiggle,&order);
704 if (wiggle<0) { wiggle=num_wiggles-1; order=wiggles[wiggle]->num_orders-1;}
705 else if (order<0 && wiggles[wiggle]->num_orders==1) order=0;
706 op=wiggles[wiggle]; nl();
707
708 if (order<0) {
709 sf(ps,SAVED_WIGGLE_WIG,maybe_s(op->num_orders),wiggle+1,1); pr();
710 if (op->num_orders>1)
711 { sf(ps,SAVED_WIGGLE_WIGS,wiggle+1,op->num_orders); pr(); }
712 sf(ps,SAVED_WIGGLE_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
713 print_old_seq(op->seq_string); nl();
714 print_saved_wiggle(wiggle);
715 } else {
716 sf(ps,SAVED_WIGGLE_WIG,"",wiggle+1,order+1); pr();
717 sf(ps,SAVED_WIGGLE_TRAIT,op->trait+1,raw.trait_name[op->trait]); pr();
718 print_old_seq(op->seq_string);
719 sf(ps,THRESHOLD_AND_SCALE,threshold,scale); pr();
720 nl(); print_wiggle_left_seq(op->data[order][0]->map);
721 last= op->data[order][0]->map->num_intervals-1;
722 print_wiggle_genetics(&op->data[order][0]->map->constraint[last]);
723 nl(); print_saved_wiggle_order(wiggle,order,threshold,scale);
724 }
725 }
726
727
728 #define WHICH_WIGGLE \
729 "%d is an ambiguous scan number.\nUse a number from %d.1-%d.%d."
730 #define PEAKS_TITLE "LOD score peaks for scan %d.%d of trait %d (%s).\n"
731 #define PEAKS_PARAMS "Peak Threshold: %-4.2lf Falloff: %-4.2lf\n"
732 #define PEAKS_NAME "Peak name %s1 has been set.\n"
733 #define PEAKS_NAMES "Peak names %s1-%s%d have been set.\n"
734 #define PEAKS_NONE "No peaks with a LOD score above the threshold were found!\n(Existing names remain unchanged).\n"
735
show_peaks()736 command show_peaks()
737 {
738 char arg[TOKLEN+1],*name;
739 int wiggle, order, last, peak;
740 real threshold, falloff;
741 WIGGLE_PEAK *peaks=NULL, *p;
742
743 name = NULL;
744 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,QCTM); /* find_wiggle_peaks calls QCTM */
745 get_arg(stoken,"",arg);
746 get_arg(rtoken,2.0,&threshold);
747 get_arg(rtoken,-2.0,&falloff);
748 nomore_args(num_args);
749
750 if (threshold<0.0 || falloff==0.0 || falloff>threshold)
751 usage_error(num_args);
752 get_wiggle_nums(arg,&wiggle,&order);
753 if (wiggle<0) { wiggle=num_wiggles-1; order=wiggles[wiggle]->num_orders-1;}
754 else if (order<0 && wiggles[wiggle]->num_orders==1) order=0;
755 if (order<0) { sf(ps,WHICH_WIGGLE,wiggle+1,wiggle+1,wiggle+1,
756 wiggles[wiggle]->num_orders-1); error(ps); }
757
758 run {
759 array(name, NAME_LEN+1, char);
760 peaks=find_wiggle_peaks(wiggle,order,threshold,-2.0,falloff,2.0,TRUE);
761 nl();
762 sf(ps,PEAKS_TITLE,wiggle+1,order+1,wiggles[wiggle]->trait+1,
763 raw.trait_name[wiggles[wiggle]->trait]); pr();
764 print_old_seq(wiggles[wiggle]->seq_string);
765 print_wiggle_left_seq(wiggles[wiggle]->data[order][0]->map);
766 last= wiggles[wiggle]->data[order][0]->map->num_intervals-1;
767 print_wiggle_genetics(
768 &wiggles[wiggle]->data[order][0]->map->constraint[last]);
769 sf(ps,PEAKS_PARAMS,threshold,falloff); pr(); nl();
770
771 for (peak=0, p=peaks; p!=NULL; p=p->next, peak++)
772 print_peak(p,peak);
773 if (peak>0) {
774 print_map_divider(); nl();
775 nstrcpy(name,raw.trait_name[wiggles[wiggle]->trait],7);
776 strcat(name,".");
777 name_peaks(peaks,name,TRUE);
778 if (peak==1) sf(ps,PEAKS_NAME,name);
779 else sf(ps,PEAKS_NAMES,name,name,peak+1);
780 } else print(PEAKS_NONE);
781
782 } when_aborting { free_wiggle_peaks(peaks); unarray(name, char) }
783 }
784
785
new_show_map()786 command new_show_map()
787 {
788 int perm,i,last,first = TRUE;
789 real inc,dist;
790
791 qtl_ready(ANY_DATA,WIGSEQ,TRAIT,QCTM);
792
793 inc = 999.0; last= -10; /* just one per interval */
794
795 print("\nlinkage maps:\n");
796 print(MAP_DIVIDER);
797 for_wiggle_orders(ints,map,inc,perm) {
798 if(map->left[0] != last && !first)
799 print(MAP_DIVIDER);
800 first = FALSE;
801 for(i = 0; i < map->num_intervals; i++) {
802 dist = map_length(map->left[i], map->right[i]);
803 if(dist <= .490) {
804 sf(ps," %s %5.1lf cM %4.1lf %%\n",
805 interval_str(map->left[i],map->right[i],TRUE),
806 haldane_cm(dist),dist*100.0);
807 pr();
808 }
809 else {
810 print(MAP_DIVIDER);
811 print("MAP:\n\n");
812 }
813 }
814 last = map->right[i-1];
815 }
816 print(MAP_DIVIDER);
817 nl();
818 }
819
820
821 #define TEST_TITLE "Test genetics results for trait %d (%s).\n"
822 #define TEST_PARAMS "Scan numbers: %d.%d-%d.%d Threshold: %-4.2lf\n"
823 #define TEST_ISNT \
824 "Scan number %d.x did not use test genetics in the scanned interval."
825
show_test_wiggle()826 command show_test_wiggle()
827 {
828 char arg[TOKLEN+1];
829 int wiggle, order;
830 real threshold;
831
832 qtl_ready(INTERCROSS,NOSEQ,NOTRAIT,NOQCTM);
833 get_arg(stoken,"",arg);
834 get_arg(rtoken,2.0,&threshold);
835 nomore_args(num_args);
836
837 if (threshold<0.0) usage_error(2);
838 get_wiggle_nums(arg,&wiggle,&order);
839 if (wiggle<0) { wiggle=num_wiggles-1; order=wiggles[wiggle]->num_orders-1;}
840 if (!isa_test_wiggle(wiggle)) error(TEST_ISNT);
841 if (order<0 && wiggles[wiggle]->num_orders==NUM_MODELS) order=0;
842 if (order<0) { sf(ps,WHICH_WIGGLE,wiggle+1,wiggle+1,wiggle+1,
843 wiggles[wiggle]->num_orders-1); error(ps); }
844
845 order= (order/NUM_MODELS)*NUM_MODELS; /* first order for this test */
846 nl();
847 sf(ps,TEST_TITLE,wiggles[wiggle]->trait+1,
848 raw.trait_name[wiggles[wiggle]->trait]); pr();
849 print_old_seq(wiggles[wiggle]->seq_string);
850 print_wiggle_left_seq(wiggles[wiggle]->data[order][0]->map);
851 sf(ps,TEST_PARAMS,wiggle+1,order+1,wiggle+1,order+NUM_MODELS,
852 threshold); pr(); nl();
853 print_test_wiggle_order(wiggle,order,threshold);
854
855 }
856
857
858 #define TRAIT_BADNAME "'%s' is not a valid name for a new trait.\n%s"
859 #define TRAITS_FULL "You have reached maximum number of traits.\nUse 'forget tr"
860 #define TRAIT_MADE "New trait number %d (%s) had been added to the data set.\n"
make_trait()861 command make_trait()
862 {
863 int i,j,k,adjusted_array_size;
864 real *normal_array=NULL;
865 NORMAL_TEST *normal_check=NULL;
866 char *eqn, *name_to_check, *name;
867 char *error_message=get_temp_string();
868 int trait_redone, trait_index;
869 EQUATION **postfixed;
870
871 /* I would like to use despace(uncrunched_args), but split_arglist
872 doesn't like that. Hmmmm... */
873 name_to_check = get_temp_string();
874 array(postfixed, MAX_EQN_SIZE, EQUATION*);
875 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
876 if (nullstr(args)) usage_error(0);
877 name=args;
878 if (!split_arglist(&eqn,'=')) usage_error(1);
879 if (!nstoken(&name,sREQUIRED,name_to_check,NAME_LEN)) usage_error(2);
880
881 if (!valid_new_trait_name(name_to_check,errmsg))
882 { sf(ps,TRAIT_BADNAME,name_to_check,errmsg); error(ps); }
883
884 trait_redone = FALSE;
885 for (i=0;i<raw.n_traits;i++) {
886 if (nullstr(raw.trait_name[i])) {
887 trait_index = i;
888 trait_redone = TRUE;
889 break;
890 }
891 }
892 if (!trait_redone) {
893 if (raw.n_traits==raw.max_traits) error(TRAITS_FULL);
894 else trait_index = raw.n_traits;
895 }
896
897 run {
898 /* The variable value and lookup stuff should be changed to read
899 right out of the raw struct. */
900 eqn_init();
901 table_size = raw.n_traits; /* check bounds? */
902 for (k=0; k<raw.n_traits; k++)
903 nstrcpy(variable_table[k],raw.trait_name[k],NAME_LEN);
904 postfixed = make_equation(eqn,variable_lookup);
905
906 for (i=0; i<raw.n_indivs; i++) { /* check bounds? */
907 for (k=0; k<raw.n_traits; k++) value_table[k]=raw.trait[i][k];
908 raw.trait[i][trait_index]=
909 evaluate_equation(postfixed,value_lookup);
910 }
911 strcpy(raw.trait_name[trait_index],name_to_check);
912 strcpy(raw.trait_eqn[trait_index],eqn);
913 if (!trait_redone) raw.n_traits++;
914
915 nl(); sf(ps,TRAIT_MADE,trait_index+1,name_to_check); pr(); nl();
916 array(normal_array,raw.n_indivs,real);
917 adjusted_array_size = raw.n_indivs;
918 for (i=0,j=0; i<raw.n_indivs; i++) {
919 if(raw.trait[i][trait_index] == MISSING_PHENO) {
920 adjusted_array_size--;
921 }
922 else {
923 normal_array[j] = raw.trait[i][trait_index];
924 j++;
925 }
926 }
927 single(normal_check, NORMAL_TEST);
928 normal_check = check_normalcy(normal_array,adjusted_array_size);
929 print_normal(normal_check, -1.0);
930 nl(); print_rhisto(normal_array,adjusted_array_size);
931 } on_exit {
932 if (msg==BADEQN) { /* msg should be 0 if none sent - CHECK THIS OUT */
933 print("error in equation: "); print(eqn); nl();
934 if (BADEQN_errpos<60) { space(BADEQN_errpos+19); print("^\n"); }
935 print(BADEQN_errmsg); nl();
936 }
937 unarray(normal_array,real);
938 /* free normal_check???? */
939 /* free postfixed - it wasn't alloced here... make it a local! */
940 }
941 }
942
943
predict()944 command predict()
945 {
946 char *name= get_temp_string();
947 int trait_redone, trait_index, i, perm;
948
949 qtl_ready(INTERCROSS,SEQ,TRAIT,QCTM);
950 if (nullstr(args)) usage_error(0);
951 if (!nstoken(&args,sREQUIRED,name,NAME_LEN)) usage_error(2);
952 if (!valid_new_trait_name(name,errmsg))
953 { sf(ps,TRAIT_BADNAME,name,errmsg); error(ps); }
954
955 for (i=0, trait_redone = FALSE; i<raw.n_traits; i++)
956 if (nullstr(raw.trait_name[i]))
957 { trait_index= i; trait_redone= TRUE; break; }
958 if (!trait_redone) {
959 if (raw.n_traits==raw.max_traits) error(TRAITS_FULL);
960 else trait_index= raw.n_traits;
961 }
962
963 if (num_orders!=1)
964 error("The current sequence spcifies more than one model.\n");
965 for_all_orders(ints,map,perm) {
966 get_fixed_qtl_weights(map);
967
968 for (i=0; i<raw.n_indivs; i++) {
969 raw.trait[i][trait_index]= model_prediction(map,i);
970 }
971 }
972 nl(); sf(ps,TRAIT_MADE,trait_index+1,name); pr(); nl();
973 if (!trait_redone) raw.n_traits++;
974 strcpy(raw.trait_name[trait_index],name);
975 strcpy(raw.trait_eqn[trait_index],ints_string);
976 }
977
978
dump_scan()979 command dump_scan()
980 {
981 char arg[TOKLEN+1], *errmsg=get_temp_string();
982 char *name=get_temp_string();
983 int t[8], num, i, j;
984 FILE *fp=NULL;
985 WIGGLE_OPERATION *op;
986 WIGGLE_INTERVAL **data;
987 WIGGLE_POINT *point;
988 int mark,position,order,wiggle;
989 real threshold,scale;
990
991 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
992 get_arg(stoken,"",arg);
993 get_arg(rtoken,2.0,&threshold);
994 get_arg(rtoken,0.25,&scale);
995 get_wiggle_nums(arg,&wiggle,&order);
996 if(wiggle<0) { wiggle=num_wiggles-1;order=wiggles[wiggle]->num_orders-1;}
997 else if (order<0 && wiggles[wiggle]->num_orders==1) order = 0;
998
999 sf(name,"scan.%d",wiggle);
1000 sf(ps,"dumping to '%s'\n",name); pr();
1001 fp= open_file(name,WRITE);
1002 op= wiggles[wiggle];
1003 data = op->data[order];
1004 position = 0;
1005 mark= 0;
1006 for (i=0;i<op->num_wiggled_intervals;i++) {
1007 if (i>0 && !data[i]->contig) { mark=0;
1008 position=(int)(position/1000)+1;
1009 position=(position*1000); }
1010
1011
1012 for (j=0;j<data[i]->num_points; j++) {
1013 point=data[i]->point[j];
1014 fprintf(fp,"%lf\t %6.3lf\t %7.3lf\n",position+mark*threshold,
1015 point->qtl_weight,point->lod_score);
1016 mark++;
1017 }
1018 }
1019 close_file(fp);
1020 }
1021
dump_traits()1022 command dump_traits()
1023 {
1024 char arg[TOKLEN+1], *errmsg=get_temp_string();
1025 int t[8], num, i, j;
1026 FILE *fp=NULL;
1027
1028 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1029
1030 num=0;
1031 do {
1032 if (num==8) error("too many traits");
1033 if (!itoken(&args,iREQUIRED,&t[num]) ||
1034 !valid_trait_num(--(t[num++]))) error("bad trait");
1035 } while (!nullstr(args));
1036
1037 run {
1038 print("dumping to 'traits.dump'\n");
1039 fp= open_file("traits.dump",WRITE);
1040 fprintf(fp,"*indiv\t");
1041 for (j=0; j<num; j++) fprintf(fp,"%7s\t",raw.trait_name[t[j]]);
1042 fprintf(fp,"\n");
1043
1044 for (i=0; i<raw.n_indivs; i++) {
1045 fprintf(fp,"%-4d\t",i+1);
1046 for (j=0; j<num; j++) {
1047 if (raw.trait[i][t[j]] != MISSING_PHENO)
1048 fprintf(fp,"%7.4lf\t",raw.trait[i][t[j]]);
1049 else fprintf(fp,"\t");
1050 }
1051 fprintf(fp,"\n");
1052 }
1053
1054 } on_exit {
1055 close_file(fp);
1056 }
1057 }
1058
1059
dump_genome()1060 command dump_genome()
1061 {
1062 FILE *fp=NULL;
1063 int i,j,k,l;
1064 int tot_aa=0,tot_bb=0,tot_ab=0,total_num=0;
1065 run {
1066 print("dumping to 'genome.dump'\n");
1067 fp = open_file("genome.dump",WRITE);
1068 fprintf(fp,"*loci\t \%AA\t \%AB \%BB\t n_indivs");
1069 for(k=0;k<raw.n_loci;k++) {
1070 fprintf(fp,"\n");
1071 for(l=0;l<raw.n_indivs;l++) {
1072 if (raw.locus[k][l]=='A') { tot_aa++; total_num++; }
1073 else if (raw.locus[k][l]=='B') {tot_bb++; total_num++; }
1074 else if (raw.locus[k][l]=='H') {tot_ab++; total_num++; }
1075 }
1076 fprintf(fp,"%s\t %lf\t %lf\t %lf\t %d",raw.locus_name[k],
1077 100*tot_aa/total_num,100*tot_bb/total_num,100*tot_ab/total_num,
1078 total_num);
1079 tot_aa=tot_bb=tot_ab=total_num=0;
1080 }
1081 fprintf(fp,"\n");
1082 fprintf(fp,"\n");
1083
1084 } on_exit {
1085 close_file(fp);
1086
1087 }
1088 }
1089
list_traits()1090 command list_traits()
1091 {
1092 int printed, i;
1093
1094 printed = 0;
1095 print(BIG_DIVIDER); print("TRAITS:\n\n");
1096 for(i=0; i< raw.n_traits; i++) {
1097 if (nullstr(raw.trait_name[i])) sf(ps, "%3d <deleted> ",i+1);
1098 else sf(ps,"%3d %-10s %s ",i+1,raw.trait_name[i],raw.trait_eqn[i]);
1099 pr();nl();
1100 }
1101 print(BIG_DIVIDER);
1102 }
1103
1104
1105 #define SHOW_WHAT "No trait specified.\nEither supply a trait name or number or use the 'trait' command."
1106
show_trait()1107 command show_trait()
1108 {
1109 real *normal_array=NULL;
1110 NORMAL_TEST *normal_check=NULL;
1111 int i, t,j, adjusted_array_size;
1112 char arg[TOKLEN+1], *errmsg=get_temp_string();
1113
1114 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1115 get_one_arg(stoken,"",arg);
1116 if (!nullstr(arg)) { if (!valid_trait_str(arg,&t,errmsg)) error(errmsg); }
1117 else { if (trait<0) error(SHOW_WHAT); else t=trait; }
1118
1119 run {
1120 sf(ps,"\nTrait %d (%s)",t+1,raw.trait_name[t]); pr();
1121 if (!nullstr(raw.trait_eqn[t]))
1122 { sf(ps,"='%s'",raw.trait_eqn[t]); pr(); }
1123 print(":\n\n");
1124 array(normal_array,raw.n_indivs,real);
1125 adjusted_array_size = raw.n_indivs;
1126 for (i=0,j=0; i<raw.n_indivs; i++) {
1127 if(raw.trait[i][t] == MISSING_PHENO) {
1128 adjusted_array_size--;
1129 }
1130 else {
1131 normal_array[j] = raw.trait[i][t];
1132 j++;
1133 }
1134 }
1135 single(normal_check, NORMAL_TEST);
1136 normal_check = check_normalcy(normal_array,adjusted_array_size);
1137 print_normal(normal_check, -1.0);
1138 nl(); print_rhisto(normal_array,adjusted_array_size);
1139 } on_exit {
1140 unarray(normal_array,real);
1141 /* free normal_check???? */
1142 }
1143 }
1144
1145
1146 #define TRAIT_NOTMADE "Can't delete specified trait.\nTrait %d (%s) was made using the 'make trait' command."
1147 #define REALLY_KILL "Trait %d (%s)"
1148
forget_trait()1149 command forget_trait()
1150 {
1151 char trait_name[TOKLEN+1], c;
1152 int t;
1153
1154 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1155 get_one_arg(stoken,sREQUIRED,trait_name);
1156 if (!valid_trait_str(trait_name,&t,errmsg)) error(errmsg);
1157 if (raw.n_traits==1) error("There is only one trait in the data set!\n.");
1158
1159 sf(ps,"Deleting trait %d (%s)",t+1,raw.trait_name[t]); pr();
1160 if (!nullstr(raw.trait_eqn[t])) { sf(ps,"='%s'",raw.trait_eqn[t]); pr(); }
1161 print("\n");
1162 getln("Are you sure you want to delete it? [no] ");
1163 if (!parse_char(&ln,"y",TRUE,&c)) return;
1164
1165 raw.trait_name[t][0]='\0'; raw.trait_eqn[t][0]='\0';
1166 if (t==raw.n_traits-1) raw.n_traits--;
1167 print("Trait Deleted.\n");
1168 if (t==trait) { trait= -1; print("The current trait is now not set!\n"); }
1169 /* need to blow away saved wiggles for this trait */
1170 }
1171 #define COMP_LIST_NAMES "%3d.%s %-10s %s\n"
1172 #define COMP_LIST_NUMS "%3d.%s %3d %s\n"
1173 #define COMP_LIST_TITLE " NUM TRAIT%s SEQUENCE\n"
1174 #define WIG_LIST_TITLE " NUM TRAIT%s SEQUENCE\n"
1175 #define WIG_LIST_NAMES "%3d.%s %-10s %s\n"
1176 #define WIG_LIST_NUMS "%3d.%s %3d %s\n"
1177
forget_compare()1178 command forget_compare()
1179 {
1180 int i,comp_number;
1181 char c;
1182
1183 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1184 get_one_arg(itoken,iREQUIRED,&comp_number);
1185 if (comp_number > num_compares) error(errmsg);
1186 sf (ps,"Deleting compare %d \n",comp_number);pr();
1187 i=comp_number-1;
1188 sf(ps,COMP_LIST_TITLE,(print_names ? " ":"")); pr();
1189 if(print_names)
1190 sf(ps,COMP_LIST_NAMES,i+1,(compares[i]->num_contigs==1
1191 ? "1" : "x"),
1192 raw.trait_name[compares[i]->trait],
1193 compares[i]->seq_string);
1194 else
1195 sf(ps,COMP_LIST_NUMS,i+1,(compares[i]->num_contigs==1
1196 ? "1" : "x"),
1197 compares[i]->trait+1,compares[i]->seq_string);
1198 pr();
1199 print("\n");
1200 getln("Are you sure you want to delete it? [no] ");
1201 if (!parse_char(&ln,"y",TRUE,&c)) return;
1202 bash_compare_struct(i);
1203 }
1204
1205
forget_all_compares()1206 command forget_all_compares()
1207 {
1208 int i;
1209 char c;
1210
1211 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1212 sf (ps, "\nDeleting ALL compares!\n");pr();
1213 getln("\nAre you sure you want to delete them? [no] ");
1214 if (!parse_char(&ln,"y",TRUE,&c)) return;
1215 for (i=0;i<num_compares;i++) {
1216 if (compares[i]->data!=NULL)
1217 bash_compare_struct(i);
1218 }
1219 num_compares = 0;
1220 }
1221
1222
forget_scan()1223 command forget_scan()
1224 {
1225 char c;
1226 int scan_number, t;
1227 WIGGLE_OPERATION *op;
1228
1229 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1230 get_one_arg(itoken,iREQUIRED,&scan_number);
1231 if (scan_number > num_wiggles) error(errmsg);
1232 sf(ps,"Deleting scan %d \n",scan_number);pr();
1233 t=scan_number-1;
1234 sf(ps,WIG_LIST_TITLE,(print_names ? " ":""));pr();
1235 op=wiggles[t];
1236 if (print_names)
1237 sf(ps,WIG_LIST_NAMES,t+1,(op->num_orders==1 ? "1":"x"),
1238 raw.trait_name[op->trait],op->seq_string);
1239 else sf(ps,WIG_LIST_NUMS,t+1,(op->num_orders==1 ? "1":"x"),
1240 op->trait+1,op->seq_string);
1241 pr();
1242 print("\n");
1243 getln("Are you sure you want to delete it? [no] ");
1244 if (!parse_char(&ln,"y",TRUE,&c)) return;
1245 bash_wiggle_struct(t);
1246 print("Scan results deleted.\n");
1247 }
1248
forget_all_scans()1249 command forget_all_scans()
1250 {
1251 char c;
1252 int i;
1253
1254 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1255 sf (ps, "Deleting ALL saved scan results!\n");pr();
1256 getln("\nAre you sure you want to delete them? [no] ");
1257 if (!parse_char(&ln,"y",TRUE,&c)) return;
1258 for (i=0;i<num_wiggles;i++) {
1259 if (wiggles[i]->data != NULL)
1260 bash_wiggle_struct(i);
1261 }
1262 num_wiggles = 0;
1263 print("All scan results deleted.\n");
1264 }
1265
1266
let()1267 command let()
1268 {
1269 char *name,*seqnce,*err;
1270
1271 name = args;
1272 split_arglist(&seqnce,'=');
1273
1274 if(!name_sequence(name,seqnce,&err))
1275 error(err);
1276 }
1277
names()1278 command names()
1279 {
1280 /* This is a KLUDGE for now until we write a macro which ports. */
1281 for (Te=context[active_context]->named_sequences->list; Te!=NULL;
1282 Te=Te->next) {
1283 sf(ps,"%-12s= %s\n",Te->id.name,Te->string);
1284 print(ps);
1285 }
1286 }
1287
forget()1288 command forget()
1289 {
1290 char *name;
1291 int fail;
1292
1293 name = get_temp_string();
1294 if(stoken(&args,sREQUIRED,name)) {
1295 if(!delete_named_entry(name,context[active_context]->named_sequences,&fail)) {
1296 if(fail == NAME_DOESNT_MATCH)
1297 sf(ps,"%s is not a defined name",name);
1298 else
1299 sf(ps,"%s is an ambiguous name",name);
1300 error(ps);
1301 }
1302 }
1303 }
1304
1305
1306 #define BADEDITEDSEQ \
1307 "An illegal sequence was specified.\nThe interval list remains '%s'.\n"
1308
sequence_editor()1309 command sequence_editor()
1310 {
1311 char name[TOKLEN+1], prompt[TOKLEN+1], *seq, *err, *foo, *new_seq;
1312 int errpos;
1313
1314 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1315 new_seq=NULL;
1316
1317 run {
1318 /* edit line will send an error if it can't be run now */
1319 array(new_seq, SEQ_LEN+1, char); new_seq[0]='\0';
1320
1321 strcpy(prompt,"sequence= ");
1322 edit_line(prompt,new_seq,SEQ_LEN,ints_string);
1323 if (!set_qtl_sequence(new_seq,errmsg,&errpos)) {
1324 print("Error in sequence '"); print(new_seq); print("'\n ");
1325 space((errpos+20) % LINE); print("^\n"); print(errmsg); nl();
1326 send(BADSEQ);
1327 }
1328 print("ok\n");
1329 unarray(new_seq,char);
1330 } except_when(BADSEQ) { unarray(new_seq,char); }
1331 }
1332
1333
1334 #define LOADED_DATA \
1335 "data files '%s' and '%s' are %sloaded.\n(%d %s progeny, %d loci, %d trait%s)\n"
1336 #define ALREADY_LOADED_DATA \
1337 "data are already loaded\nshut down and restart the program to load new data"
1338 #define NO_LOADED_DATA \
1339 "no data are loaded\ntype \"load data <filename>\" to load data\n"
1340
load_data()1341 command load_data()
1342 {
1343 bool prev_dat, unloaded;
1344 FILE *fpa=NULL, *fpb=NULL, *fpc=NULL;
1345 int num_of_file=0;
1346 char *dfile= get_temp_string(), *tfile= get_temp_string(), *mfile= get_temp_string();
1347
1348 if (nullstr(args)) {
1349 if (!data_loaded()) print(NO_LOADED_DATA); else {
1350 strcpy(tfile,raw.file);
1351 make_filename(tfile,FORCE_EXTENSION,WRS(TRAIT_EXT));
1352 sf(ps,LOADED_DATA,raw.file,tfile,"",raw.n_indivs,
1353 (raw.data_type==BACKCROSS ? "backcross":"intercross"),
1354 raw.n_loci,raw.n_traits,maybe_s(raw.n_traits)); pr();
1355 }
1356
1357 } else { /* !nullstr */
1358
1359 /* for now we prevent reloading */
1360 if ((prev_dat = data_loaded())) error(ALREADY_LOADED_DATA);
1361
1362 unloaded=FALSE;
1363 nstoken(&args,sREQUIRED,dfile,PATH_LENGTH); tfile[0]='\0';
1364 nomore_args(num_args);
1365 run {
1366 if (!make_filename(dfile,FORCE_EXTENSION,WRS(DATA_EXT))) send(CANTOPEN);
1367 fpa= open_file(dfile,READ);
1368 fgetln_(fpa);
1369
1370 if (streq(ln,"prepared data f2 backcross")) {
1371 strcpy(geno_chars,default_backcross_chars);
1372 raw.data_type=BACKCROSS;
1373
1374 } else if (streq(ln,"prepared data f2 intercross")) {
1375 strcpy(geno_chars,default_intercross_chars);
1376 raw.data_type=INTERCROSS; raw.f3=FALSE;
1377
1378 } else if (streq(ln,"prepared data f3")) {
1379 strcpy(geno_chars,default_intercross_chars);
1380 raw.data_type=INTERCROSS; raw.f3=FALSE; /* KLUDGE!!! */
1381
1382 } else {
1383 error("unrecognized header line in data file");
1384 }
1385
1386 strcpy(tfile,dfile);
1387 if (!make_filename(tfile,FORCE_EXTENSION,WRS(TRAIT_EXT)))
1388 send(CANTOPEN);
1389 fpb = open_file(tfile,READ);
1390 strcpy(mfile,dfile);
1391 if (!make_filename(mfile,FORCE_EXTENSION,WRS(MAPS_EXT)))
1392 send(CANTOPEN);
1393 fpc = open_file(mfile,READ);
1394
1395 read_data(fpa,fpb,fpc,dfile,num_of_file);
1396
1397 crunch_data();
1398 allocate_qtl_struct(raw.max_traits*2,raw.max_traits*2);
1399
1400 sf(ps,LOADED_DATA,dfile,tfile,"",raw.n_indivs,
1401 (raw.data_type==BACKCROSS ? "backcross":"intercross"),
1402 raw.n_loci-dum_loc,raw.n_traits,maybe_s(raw.n_traits)); pr();
1403 if (raw.n_traits==1) trait=0; else trait= NOTRAIT;
1404 update_top();
1405 if(!altered_chroms)
1406 load_qtl_files();
1407
1408 } on_exit {
1409 if (msg==CANTOPEN) {
1410 sf(ps,"error: unable to open data file\n"); pr();
1411 } else if (msg==BADDATA) {
1412 sf(ps,"error: unable to load data from file\nline %d:",
1413 temp,BADDATA_line_num); pr();
1414 print(BADDATA_error); nl();
1415 strcpy(raw.file,"");
1416 }
1417 if (msg!=CANTOPEN && msg!=0) {
1418 print("data file not loaded\n");
1419 strcpy(raw.file,"");
1420 close_file(fpa); close_file(fpb); close_file(fpc);
1421 relay_messages;
1422 }
1423 }
1424 }
1425 }
1426
1427 /********** NO MORE PREP IN QTL
1428
1429
1430 command prep_data()
1431 {
1432 char *in_name,*out_name,*str;
1433 FILE *fp;
1434 int length;
1435
1436 in_name = NULL; out_name = NULL;
1437
1438 run {
1439 str = get_temp_string();
1440 if(!stoken(&args,sREQUIRED,str))
1441 input("Input file to be prepared: ",str,25);
1442 in_name = str;
1443
1444 length = len(str);
1445 if(streq(&str[length-5],DATA_EXT))
1446 error("raw file cannot have '.data' extension");
1447 if(streq(&str[length-7],TRAIT_EXT))
1448 error("raw file cannot have '.trait' extension");
1449
1450 out_name = get_temp_string();
1451 strcpy(out_name,in_name);
1452 make_filename(out_name,FORCE_EXTENSION,WRS(DATA_EXT));
1453
1454 fp = open_file(in_name,READ);
1455 getdataln(fp); crunch(ln);
1456 close_file(fp);
1457 if(streq(ln,"data type f2 intercross") ||
1458 streq(ln,"data type f3 intercross") ||
1459 streq(ln,"data type f2 backcross") ||
1460 streq(ln,"data type ri sib") ||
1461 streq(ln,"data type ri self") ||
1462 streq(ln,"data type f2") ||
1463 streq(ln,"data type backcross") ||
1464 streq(ln,"data type quant intercross") ||
1465 streq(ln,"data type quant backcross") ) {
1466 f2_prep(in_name,out_name);
1467 }
1468 else {
1469 sf(ps,"'%s' is not an allowable data type",ln);
1470 error(ps);
1471 }
1472 } except_when(CANTOPEN) {
1473 sf(ps,"Can't open file '%s'.",in_name);
1474 error(ps);
1475
1476 }
1477 }
1478
1479 **********/
1480
save_status()1481 command save_status()
1482 {
1483 char *name, *name2, *name3;
1484 FILE *fp;
1485
1486 if(nullstr(raw.file)) {
1487 print("No data have been loaded, none can be saved.\n"); return; }
1488 name = mkstrcpy(raw.file);
1489 name3 = mkstrcpy(raw.file);
1490 name2 = "dummy.qtls";
1491 make_filename(name, FORCE_EXTENSION, WRS(QTL_EXT));
1492
1493 /* The following accounts for the case where there is no file 'name'
1494 (i.e. it makes one) */
1495 run {
1496 fp = open_file(name,READ);
1497 close_file(fp);
1498 } except_when(CANTOPEN) {
1499 fp = open_file(name,WRITE);
1500 close_file(fp);
1501 }
1502
1503 make_filename(name3, FORCE_EXTENSION, WRS(QTL_OLD));
1504 run {
1505 fp = open_file(name2, WRITE);
1506 sf(ps,"Now saving %s...\n",name); pr();
1507 save_qtl_files(fp);
1508 close_file(fp);
1509 if (rename_file(name,name3)) rename_file(name2,name);
1510 } except {
1511 when CANTOPEN:
1512 sf(ps,"Can't open %s.\n",name); pr();
1513 default:
1514 rename_file(name3,name);
1515 if(msg == INTERRUPT) send(INTERRUPT);
1516 print("\nAn error occurred while the data files were being rewritten.\n");
1517 print("You may have to turn the 'auto save' option 'off' before quitting.\n");
1518 if(redirecting_input) send(QUIT);
1519
1520 }
1521
1522 name2="dummy.traits";
1523 make_filename(name, FORCE_EXTENSION, WRS(TRAIT_EXT));
1524 make_filename(name3, FORCE_EXTENSION, WRS(TRAIT_OLD));
1525 run {
1526 fp = open_file(name2, WRITE);
1527 sf(ps,"Now saving %s...\n",name); pr();
1528 save_traitfile(fp);
1529 close_file(fp);
1530 if (rename_file(name,name3)) rename_file(name2,name);
1531 } except {
1532 when CANTOPEN:
1533 sf(ps,"Can't open fle '%s'.\n",name); pr(); /* fall through */
1534 default:
1535 rename_file(name3,name);
1536 if(msg == INTERRUPT) send(INTERRUPT);
1537 print("\nAn error occurred while the data files were being rewritten.\n");
1538 print("You may have to turn the 'auto save' option 'off' before quitting.\n");
1539 if (redirecting_input) send(QUIT);
1540 }
1541 }
1542
1543
save_on_exit(do_it_now)1544 bool save_on_exit(do_it_now)
1545 bool do_it_now;
1546 {
1547 if(!do_it_now) return(auto_save && data_loaded());
1548 if(auto_save && data_loaded()) save_status();
1549 return(TRUE);
1550 }
1551
save_qtl_files(fp)1552 void save_qtl_files(fp)
1553 FILE *fp;
1554 {
1555 int i;
1556
1557 fprintf(fp,"%d\n",raw.filenumber);
1558 fprintf(fp,"%d %d\n",num_wiggles,first_wiggle);
1559 for(i = 0; i < num_wiggles; i++) {
1560 save_wiggle(fp,i);
1561 }
1562
1563 fprint(fp,WRS("#Compares\n"));
1564
1565 fprintf(fp,"%d %d\n",num_compares,first_compare);
1566 for(i = 0; i < num_compares; i++) {
1567 save_compare(fp,i);
1568 }
1569 }
1570
1571 #define QTL_LD_ERROR1 \
1572 "An error ocurred while loading the qtl map data in file '%s'.\n"
1573 #define QTL_LD_ERROR2 \
1574 "In some extreme cases, this can prevent correct operation of this program.\n"
1575 #define QTL_LD_ERROR3 \
1576 "If problems occur, delete this file and try reloading the data set.\n"
1577 #define QTL_LD_OLD1 \
1578 "You have prepared a new version of your data files since the last qtl map\n"
1579 #define QTL_LD_OLD2 \
1580 "data were saved in the '.qtls' file. These data will not be loaded.\n"
1581
load_qtl_files()1582 void load_qtl_files()
1583 {
1584 char name[PATH_LENGTH+1];
1585 FILE *fp=NULL;
1586 int i, n_wigs, n_comps, filenum;
1587
1588 run {
1589 strcpy(name,raw.file);
1590 make_filename(name,FORCE_EXTENSION,WRS(QTL_EXT));
1591 fp= open_file(name,READ);
1592
1593 if(fscanf(fp,"%d\n",&filenum) != 1) send(IOERROR); /* KLUDGE */
1594 else if(filenum != raw.filenumber) {
1595 print(QTL_LD_OLD1); print(QTL_LD_OLD2);
1596 }
1597 else if(fscanf(fp,"%d %d\n",&n_wigs,&first_wiggle)!= 2) send(IOERROR);
1598
1599 else {
1600 for (i = 0; i < n_wigs; i++) {
1601 load_wiggle(fp);
1602 }
1603
1604 fgetdataln(fp,NULL);
1605
1606 if (sscanf(ln,"%d %d\n",&n_comps,&first_compare) != 2)
1607 send(IOERROR);
1608
1609 for(i = 0; i < n_comps; i++) {
1610 load_compare(fp);
1611 }
1612
1613 sf(ps,"QTL map data in file '%s' have been loaded.\n",name);
1614 pr();
1615 }
1616 close_file(fp);
1617
1618 } except {
1619 when CANTOPEN:
1620 print("Unable to load any saved QTL map data.\n");
1621 msg=0; break; /* No .qtls file - WHY IS THIS MSG=0 NEEDED??? */
1622 default:
1623 sf(ps,QTL_LD_ERROR1,name); pr();
1624 print(QTL_LD_ERROR2); print(QTL_LD_ERROR3);
1625 close_file(fp);
1626 }
1627 }
1628
1629
1630
tester()1631 command tester()
1632 {
1633 real theta, pos, f2_sum, f3_sum, a, b, c, x, y, z, left_rf, right_rf;
1634 int geno, qtl, left, right;
1635
1636 getln("Theta, L_Pos, R_pos: ");
1637 sscanf(ln,"%lf %lf %lf",&theta,&left_rf,&right_rf);
1638 /* left_rf= rmaxf(MIN_REC_FRAC,min(pos,MAX_FRAC_OF_RF*theta));
1639 right_rf= (theta - left_rf)/(1 - 2*left_rf); */
1640
1641 sf(ps,"left_rf=%lf right_rf=%lf\n",left_rf,right_rf); pr();
1642 for (left=0; left<4; left++)
1643 for (right=0; right<4; right++) {
1644 f2_sum=f3_sum= 0.0;
1645 x=y=z= 0.0;
1646 raw.f3=TRUE; c=transition_prob(INTERCROSS,left,right,theta);
1647 /* raw.f3=FALSE; z=transition_prob(INTERCROSS,left,right,theta); */
1648 sf(ps,"left=%d right=%d F2-prob=%lf F3-prob=%lf\n",left,right,z,c); pr();
1649 for (qtl=0; qtl<4; qtl++) {
1650 raw.f3=TRUE; f3_sum+=
1651 ((a=transition_prob(INTERCROSS,left,qtl,left_rf)) *
1652 (b=transition_prob(INTERCROSS,qtl,right,right_rf))) /c;
1653 /* raw.f3=FALSE; f2_sum+=
1654 ((x=transition_prob(INTERCROSS,left,qtl,left_rf)) *
1655 (y=transition_prob(INTERCROSS,qtl,right,right_rf))) /z; */
1656 sf(ps,"qtl= %d F2: L=%lf R=%lf F3: L=%lf R=%lf\n",qtl,x,y,a,b); pr();
1657 }
1658 sf(ps,"==== F2 sum= %lf F3 sum=%lf\n",f2_sum,f3_sum); pr();
1659 }
1660 }
1661
1662
tweak_weight()1663 command tweak_weight()
1664 {
1665 real start,end,step,weight;
1666 int perm;
1667
1668 qtl_ready(ANY_DATA,SEQ,TRAIT,QCTM);
1669 get_arg(rtoken,0.0,&start);
1670 get_arg(rtoken,1.0,&end);
1671 get_arg(rtoken,0.1,&step);
1672 nomore_args(num_args);
1673
1674 run {
1675 fix_weight_kludge = TRUE;
1676 for_all_orders(ints,map,perm) {
1677 nl(); print_wiggle_title();
1678 print_wiggle_interval(map);
1679 for (weight=start; weight<=end; weight+=step) {
1680 map->qtl_weight[0]= weight;
1681 make_qtl_map(map);
1682 print_wiggle_map(map,2.0,0.25);
1683 }
1684 print_wiggle_interval(NULL);
1685 }
1686 } on_exit { fix_weight_kludge = FALSE; }
1687 fix_weight_kludge = FALSE;
1688 }
1689
1690
draw_wiggle()1691 command draw_wiggle()
1692 {
1693 char arg[TOKLEN+1];
1694 int wiggle, order, last;
1695 real threshold, scale;
1696 WIGGLE_OPERATION *op;
1697
1698 qtl_ready(ANY_DATA,NOSEQ,NOTRAIT,NOQCTM);
1699 get_arg(stoken,"",arg);
1700 get_arg(rtoken,2.0,&threshold);
1701 get_arg(rtoken,0.25,&scale);
1702 nomore_args(num_args);
1703
1704 if (threshold<0.0 || scale<=0.0) usage_error(2);
1705 get_wiggle_nums(arg,&wiggle,&order);
1706 if (wiggle<0) { wiggle=num_wiggles-1; order=wiggles[wiggle]->num_orders-1;}
1707 else if (order<0 && wiggles[wiggle]->num_orders==1) order=0;
1708 op=wiggles[wiggle]; nl();
1709
1710 if (order<0) {
1711 print_ps_multi_wiggle(wiggle, threshold);
1712 } else {
1713 print_ps_wiggle_order(wiggle, order, threshold);
1714 }
1715 }
1716
1717
1718