1 /******************************************************************************
2
3 #### # # # #### #### # ###### ####
4 # # # # # # # # # # # # #
5 # # # # # # # # ##### #
6 # # # # ## # # # ### # ### # # ### #
7 # # ## ## # # # # # # # ### # #
8 ### # # # # #### #### ###### ###### ### ####
9
10 ******************************************************************************/
11
12 /* qwiggle.c - code for global wiggle/compare data storage, save & load */
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 /* globals */
22 /* WIGGLE_INTERVAL ****qtls; */
23 WIGGLE_OPERATION **wiggles;
24 int num_wiggles, max_wiggles, first_wiggle;
25 COMPARE_OPERATION **compares;
26 int num_compares, max_compares, first_compare;
27
28 /* local to this file */
29 bool step();
30 WIGGLE_PEAK *peak_finder();
31 void save_wiggle(),save_interval(),save_qtl_map(),save_compare();
32 void load_wiggle(),load_interval(),load_compare();
33 QTL_MAP *load_qtl_map();
34
35 /* functions */
36
wiggle_init()37 void wiggle_init()
38 {
39 /* qtls=NULL; */
40 wiggles=NULL;
41 compares=NULL;
42 num_wiggles= num_compares= 0;
43 first_wiggle= first_compare= 0;
44 max_wiggles= max_compares= 0; /* WHAT DEFAULT TO USE? */
45 }
46
47
48
allocate_qtl_struct(n_wiggles,n_compares)49 void allocate_qtl_struct(n_wiggles,n_compares) /* allocate globals */
50 int n_wiggles, n_compares;
51 /* use globals raw.max_traits, raw.data_type, and raw.n_loci as params */
52 {
53 int i, j, k, n_models, n_intervals;
54
55 if (raw.data_type==NO_DATA) send(CRASH);
56 wiggles=NULL; compares=NULL; /* qtls=NULL; */
57
58 run {
59 parray(wiggles,n_wiggles,WIGGLE_OPERATION);
60 max_wiggles=n_wiggles; num_wiggles=0;
61 for (i=0; i<n_wiggles; i++) wiggles[i]->data=NULL;
62
63 parray(compares,n_compares,COMPARE_OPERATION);
64 max_compares=n_compares; num_compares=0;
65 for (i=0; i<n_compares; i++) compares[i]->data=NULL;
66
67 n_models= (raw.data_type==INTERCROSS ? NUM_MODELS:1);
68 n_intervals= raw.n_loci-1;
69 /* matrix(qtls,raw.n_traits,n_models,WIGGLE_INTERVAL**);
70 for (i=0; i<raw.n_traits; i++) for (j=0; j<n_models; j++) {
71 array(qtls[i][j],n_intervals,WIGGLE_INTERVAL*);
72 for (k=0; k<n_intervals; k++) qtls[i][j][k]=NULL;
73 } */
74
75 } except_when(NOMEMORY) {
76 unparray(wiggles,n_wiggles,WIGGLE_OPERATION);
77 unparray(compares,n_compares,compare_operation);
78 /* if (qtls!=NULL)
79 for (i=0; i<raw.n_traits; i++) for (j=0; j<n_models; j++)
80 unparray(qtls[i][j],n_intervals,WIGGLE_INTERVAL);
81 unmatrix(qtls,raw.n_traits,WIGGLE_INTERVAL**); */
82 wiggles=NULL; compares=NULL; /* qtls=NULL; */
83 num_wiggles=max_wiggles=num_compares=max_compares=0;
84 relay;
85 }
86 }
87
88
89 /***** Stuff for the wiggle struct to save wiggle output... *****/
90
allocate_wiggle_struct(trait,seq,seq_str,n_intervals,n_orders,n_wiggled)91 int allocate_wiggle_struct(trait,seq,seq_str,n_intervals,n_orders,n_wiggled)
92 int trait;
93 QTL_SEQUENCE *seq;
94 char *seq_str;
95 int n_intervals, n_orders, n_wiggled;
96 /* return wiggles struct entry#, or -1 if error */
97 {
98 int n, i, j;
99
100 /* KLUDGE - How to recycle these things? */
101 if (num_wiggles==max_wiggles) return(-1);
102
103 n= num_wiggles++;
104 wiggles[n]->data=NULL;
105 wiggles[n]->seq_string=NULL;
106 run {
107 if (trait == -1) {
108 wiggles[n]->trait= -1;
109 wiggles[n]->seq = NULL;
110 wiggles[n]->num_intervals = -1;
111 wiggles[n]->num_orders = -1;
112 wiggles[n]->num_wiggled_intervals = -1;
113 wiggles[n]->order = 0;
114 wiggles[n]->wiggle_interval = -1;
115 wiggles[n]->seq_string = NULL;
116 }
117 else {
118 wiggles[n]->trait=trait;
119 wiggles[n]->seq= seq; seq->dont_free=TRUE;
120 wiggles[n]->num_intervals=n_intervals;
121 wiggles[n]->num_orders= n_orders;
122 wiggles[n]->num_wiggled_intervals= n_wiggled;
123 wiggles[n]->order= -1; /* For store_wiggle_interval */
124 wiggles[n]->wiggle_interval= 0;
125 /**********
126 would like to expand the sequence so that changes to user-defined
127 names do not ruin any saved scan information, but this can't be
128 done here - too many problems - just DON'T edit user-def names!
129 wiggles[n]->seq_string= expand_named_entries(seq_str);
130 **********/
131 wiggles[n]->seq_string= mkstrcpy(seq_str);
132
133 matrix(wiggles[n]->data,n_orders,n_wiggled,WIGGLE_INTERVAL*);
134 for (i=0; i<n_orders; i++) for (j=0; j<n_wiggled; j++) {
135 single(wiggles[n]->data[i][j],WIGGLE_INTERVAL);
136 wiggles[n]->data[i][j]->point=NULL;
137 wiggles[n]->data[i][j]->map=NULL;
138 }
139 }
140 } when_aborting { bash_wiggle_struct(n); relay; }
141 return(n);
142 }
143
144
bash_wiggle_struct(n)145 void bash_wiggle_struct(n)
146 /* this is a KLUDGE for now - we need to recycle these better */
147 int n;
148 {
149 unmatrix(wiggles[n]->data,wiggles[n]->num_orders,WIGGLE_INTERVAL);
150 unarray(wiggles[n]->seq_string,char);
151
152 /* change this if any functions besides allocate_wiggle_struct()
153 or allocate_compare_struct() could set seq->dont_fix! */
154 wiggles[n]->seq->dont_free=FALSE;
155 wiggles[n]->data=NULL; /* This is the flag that says 'not filled in' */
156 wiggles[n]->seq_string=NULL;
157 if (n==num_wiggles-1) num_wiggles--;
158 }
159
160
store_wiggle_interval(wiggle_num,map,new_left_order,contig,cm_step)161 void store_wiggle_interval(wiggle_num,map,new_left_order,contig,cm_step)
162 int wiggle_num;
163 QTL_MAP *map;
164 bool new_left_order;
165 bool contig;
166 real cm_step;
167 {
168 WIGGLE_OPERATION *op;
169 WIGGLE_INTERVAL *data, *prev;
170 int i, j, k, n, t;
171 real cm, interval_cm;
172 bool new;
173
174 if (wiggle_num<0) return;
175
176 if (wiggle_num>=max_wiggles || wiggles[wiggle_num]==NULL ||
177 (op=wiggles[wiggle_num])->data==NULL || map==NULL) send(CRASH);
178
179 if (op->order>=0) { /* prev order exists */
180 prev=op->data[op->order][op->wiggle_interval];
181 if (prev==NULL || prev->point_num!=prev->num_points-1) send(CRASH);
182 }
183
184 if (new_left_order) {
185 if (op->order>=0) {
186 if (op->order>=op->num_orders-1) send(CRASH);
187 if (op->wiggle_interval!=op->num_wiggled_intervals-1) send(CRASH);
188 }
189 j=(op->wiggle_interval=0); i=(op->order+=1);
190 } else {
191 if (op->wiggle_interval>=op->num_wiggled_intervals-1) send(CRASH);
192 i=op->order; j=(op->wiggle_interval+=1);
193 }
194
195 k= map->num_intervals-1; /* the rightmost interval */
196 data=op->data[i][j]; if (data==NULL) send(CRASH);
197 data->point=NULL; data->map=NULL;
198 run {
199 data->contig= contig;
200 data->map= alloc_qtl_map(map->num_intervals,map->num_continuous_vars);
201 mapcpy(data->map,map);
202 data->max_like_map= FALSE;
203 data->in_qtls_struct= FALSE;
204 data->op= op;
205 data->cm_increment= cm_step;
206 /* don't change this without changing the code in wiggle_perm() */
207 interval_cm= haldane_cm(map_length(map->left[k],map->right[k]));
208 for (cm=0.0, n=0; cm<interval_cm; cm+=cm_step) n++;
209 data->num_points=n;
210 parray(data->point,data->num_points,WIGGLE_POINT);
211 data->point_num= -1; /* so that store_wiggle_point can increment */
212 data->best_point=0;
213
214 } when_aborting {
215 free_qtl_map(data->map);
216 unparray(data->point,data->num_points,WIGGLE_POINT);
217 data->point= NULL; data->map=NULL;
218 relay;
219 }
220
221 /* This wiggle interval MAY want to go into the qtls struct also... */
222 do {
223 if (map->num_intervals!=1) break;
224 t= map->trait;
225 if (raw.data_type==INTERCROSS) {
226 if ((i=map->constraint[0].interx_type)>=NUM_MODELS) break;
227 } else { i=0; if (map->constraint[0].backx_weight!=DONT_FIX) break; }
228 if ((j=map->left[0])!=map->right[0]-1) break;
229 #ifdef DELETED
230 /* It does want to go in... */
231 /* if (qtls[t][i][j]!=NULL) { but these wiggle data already exist! */
232 /* prev= qtls[t][i][j]; new=FALSE;
233 /* If the old one has points, we keep the higher resolution one */
234 if (prev->point==NULL || data->cm_increment<=prev->cm_increment) {
235 prev->in_qtls_struct=FALSE;
236 qtls[t][i][j]=data;
237 data->in_qtls_struct=TRUE;
238 new=TRUE;
239 }
240 /* If the max_like map, but not the points, was filled in
241 in the old struct we merge the ML map into the new struct. */
242 if (new && prev->max_like_map) {
243 mapcpy(data->map,prev->map);
244 data->max_like_map= TRUE;
245 }
246 /* If ONLY the max_like_map exists in the old wiggle_interval
247 struct then it should be freed. */
248 if (new && prev->point==NULL) {
249 free_qtl_map(prev->map); prev->map=NULL;
250 unsingle(prev,WIGGLE_INTERVAL);
251 }
252
253 else { /* No previous data exist. */
254 qtls[t][i][j]=data;
255 data->in_qtls_struct=TRUE;
256 }
257 #endif
258 } while(FALSE); /* not really a loop, just used do {} for convenience */
259 }
260
261
store_wiggle_point(wiggle_num,map)262 void store_wiggle_point(wiggle_num,map)
263 int wiggle_num;
264 QTL_MAP *map;
265 {
266 WIGGLE_OPERATION *op;
267 WIGGLE_INTERVAL *data;
268 int i, j;
269
270 if (wiggle_num<0) return;
271
272 if (wiggle_num>=max_wiggles || wiggles[wiggle_num]==NULL ||
273 (op=wiggles[wiggle_num])->data==NULL || map==NULL) send(CRASH);
274
275 data=op->data[op->order][op->wiggle_interval];
276 if (data==NULL || data->map==NULL) send(CRASH);
277 if (data->point_num>=data->num_points-1) send(CRASH);
278 i= (data->point_num+=1);
279 j= map->num_intervals-1; /* the rightmost interval */
280
281 data->point[i]->qtl_pos= map->qtl_pos[j];
282 data->point[i]->lod_score= map->log_like;
283 data->point[i]->var_explained= map->var_explained;
284 data->point[i]->qtl_weight= map->qtl_weight[j];
285 if (raw.data_type!=INTERCROSS) data->point[i]->qtl_dominance= 0.0;
286 else data->point[i]->qtl_dominance= map->qtl_dominance[j];
287
288 /* the "|| i == 0" was added to copy the first map into the best map
289 as an initial condition - mjd 4/3/90 */
290
291 if (data->point[data->best_point]->lod_score < map->log_like || i == 0)
292 { data->best_point=i; mapcpy(data->map,map); }
293 }
294
295
296
297
298 /***** Stuff for the COMPARE struct to save compare output... *****/
299
allocate_compare_struct(trait,seq,seq_str,n_intervals,n_orders)300 int allocate_compare_struct(trait,seq,seq_str,n_intervals,n_orders)
301 int trait;
302 QTL_SEQUENCE *seq;
303 char *seq_str;
304 int n_intervals, n_orders;
305 /* return compare struct entry#, or -1 if error */
306 {
307 int n, i, j;
308
309 /* KLUDGE - How to recycle these things? */
310 if (num_compares==max_compares) return(-1);
311
312 n= num_compares++;
313 compares[n]->data=NULL; compares[n]->seq_string=NULL;
314 run {
315 if (trait== -1){
316 compares[n]->trait= -1;
317 compares[n]->seq = NULL;
318 compares[n]->num_intervals= -1;
319 compares[n]->num_orders = -1;
320 compares[n]->num_contigs = -1;
321 compares[n]->order = 0;
322 compares[n]->seq_string = NULL;
323 compares[n]->data = NULL;
324 }
325 else {
326 compares[n]->trait=trait;
327 compares[n]->seq= seq; seq->dont_free=TRUE;
328 compares[n]->num_intervals=n_intervals;
329 compares[n]->num_orders= n_orders;
330 compares[n]->num_contigs= 0;
331 compares[n]->order= -1; /* For store_compare_interval */
332 compares[n]->seq_string= mkstrcpy(seq_str);
333 parray(compares[n]->data,n_orders,COMPARE_DATA);
334 }
335 } when_aborting { bash_compare_struct(n); relay; }
336 return(n);
337 }
338
339
bash_compare_struct(n)340 void bash_compare_struct(n)
341 /* this is a KLUDGE for now - we need to recycle these better */
342 int n;
343 {
344 free_qtl_map(compares[n]->data[0]->map);
345 compares[n]->data[0]->map=NULL;
346 unparray(compares[n]->data,compares[n]->num_orders,COMPARE_DATA);
347 unarray(compares[n]->seq_string,char);
348 /* change this if any functions besides allocate_wiggle_struct()
349 or allocate_compare_struct() could set seq->dont_fix! */
350 compares[n]->seq->dont_free=FALSE;
351 compares[n]->data=NULL; /* This is the flag that says 'not filled in' */
352 compares[n]->seq_string=NULL;
353 if (n==num_compares-1) num_compares--;
354 }
355
356
store_compare_map(compare_num,map,contig)357 void store_compare_map(compare_num,map,contig)
358 int compare_num;
359 QTL_MAP *map;
360 bool contig;
361 {
362 COMPARE_OPERATION *op;
363 COMPARE_DATA *data;
364 WIGGLE_INTERVAL *wig, *prev;
365 int i, j, k, n, t;
366
367 if (compare_num<0) return;
368 if (compare_num>=max_compares || compares[compare_num]==NULL ||
369 (op=compares[compare_num])->data==NULL || map==NULL) send(CRASH);
370
371 if (op->order<0 || !contig) op->num_contigs++;
372 i=(op->order+=1); if (i>=op->num_orders) send(CRASH);
373 data=op->data[i]; if (data==NULL) send(CRASH);
374 data->map=NULL;
375
376 run {
377 data->contig= contig;
378 data->map= alloc_qtl_map(map->num_intervals,map->num_continuous_vars);
379 mapcpy(data->map,map);
380 data->op= op;
381 } when_aborting {
382 free_qtl_map(data->map);
383 data->map=NULL;
384 relay;
385 }
386
387 /* This map MAY want to go into the qtls struct also... */
388 wig=NULL;
389 do { run {
390 if (map->qtl_pos[0]!=DONT_FIX) break;
391 if (map->num_intervals!=1) break;
392 t= map->trait;
393 if (raw.data_type==INTERCROSS) {
394 if ((i=map->constraint[0].interx_type)>=NUM_MODELS) break;
395 } else { i=0; if (map->constraint[0].backx_weight!=DONT_FIX) break; }
396 if ((j=map->left[0])!=map->right[0]-1) break;
397 #ifdef DELETED
398 /* It does want to go in... */
399 if (qtls[t][i][j]!=NULL) { /* but data for this qtl already exists! */
400 prev= qtls[t][i][j];
401 /* Regardless of whether the old one has a max-like qtl_map,
402 we give it this one. */
403 prev->max_like_map= TRUE;
404 mapcpy(prev->map,data->map);
405 } else { /* No previous data exist, so we alloc a dummy wiggle_int. */
406 wig=NULL;
407 single(wig,WIGGLE_INTERVAL);
408 wig->map=NULL; wig->point=NULL;
409 wig->map=
410 alloc_qtl_map(map->num_intervals,map->num_continuous_vars);
411 mapcpy(wig->map,data->map);
412 wig->in_qtls_struct=TRUE;
413 wig->max_like_map= TRUE;
414 qtls[t][i][j]=wig;
415 }
416 #endif
417 } when_aborting {
418 if (wig!=NULL) free_qtl_map(wig->map);
419 unsingle(wig,WIGGLE_INTERVAL);
420 relay;
421 } } while(FALSE); /* not really a loop, just used do {} for convenience */
422 }
423
424
425
step(interval,point,contig,wig,n_intervals,forwards)426 bool step(interval,point,contig,wig,n_intervals,forwards)
427 int *interval, *point;
428 bool *contig;
429 WIGGLE_INTERVAL **wig;
430 int n_intervals;
431 bool forwards;
432 {
433 if (*interval >= n_intervals) return(FALSE);
434 if (forwards) {
435 if (++*point<wig[*interval]->num_points) { *contig=TRUE; return(TRUE);}
436 else if (++*interval<n_intervals) {
437 *contig= wig[*interval]->contig;
438 *point=0; return(TRUE);
439 } else {
440 *contig = FALSE;
441 return(FALSE);
442 }
443 } else { /* backwards */
444 if (--*point>=0) { *contig=TRUE; return(TRUE); }
445 else if (--*interval>=0) {
446 if(*interval+1 == n_intervals) *contig = TRUE;
447 else *contig= wig[*interval + 1]->contig;
448 *point=wig[*interval]->num_points-1; return(TRUE);
449 } else {
450 *contig = FALSE;
451 return(FALSE);
452 }
453 }
454 }
455
456
find_wiggle_peaks(wiggle_num,left_order_num,threshold,qtl_falloff,confidence_falloff,min_peak_delta,get_peak_maps)457 WIGGLE_PEAK *find_wiggle_peaks(wiggle_num,left_order_num,threshold,qtl_falloff,
458 confidence_falloff,min_peak_delta,get_peak_maps)
459 int wiggle_num, left_order_num;
460 real threshold, qtl_falloff, confidence_falloff, min_peak_delta;
461 bool get_peak_maps;
462 {
463 int i, j, peak_i, peak_j, n_intervals;
464 bool still_falling, off_end, contig;
465 real lod, local_maxima, pos, cm;
466 real local_minima, prev_lod, starting_value;
467 WIGGLE_INTERVAL **interval;
468 WIGGLE_PEAK *first, *last, *peak;
469
470 if (wiggle_num<0) return(NULL);
471
472 if (wiggles[wiggle_num]==NULL || wiggles[wiggle_num]->data==NULL ||
473 (interval=wiggles[wiggle_num]->data[left_order_num])==NULL)
474 send(CRASH);
475
476 /* Look for a local maxima over the LOD threshold... */
477 i=0; j=0;
478 contig=FALSE; first=last=NULL;
479 n_intervals= wiggles[wiggle_num]->num_wiggled_intervals;
480
481 do {
482 if (!contig) {
483 local_maxima= -10.0; local_minima= 1000.0; lod= 1000.0;
484 peak_i=i; peak_j=j; still_falling=FALSE; off_end=FALSE;
485 starting_value = interval[i]->point[j]->lod_score;
486 }
487 prev_lod=lod;
488 lod=interval[i]->point[j]->lod_score;
489 if (!still_falling && lod>local_maxima && local_maxima>threshold &&
490 (local_maxima>=local_minima+min_peak_delta ||
491 local_minima >= starting_value + confidence_falloff)) {
492 /* Found a local maxima - now get peak and confidence interval. */
493 peak=peak_finder(&peak_i,&peak_j,&off_end,get_peak_maps,
494 interval,n_intervals,threshold,
495 confidence_falloff);
496 /* Store the data for this peak */
497 if (first==NULL) first=last=peak;
498 else { last->next=peak; last=peak; }
499 /* peak_finder sets the counters to be one step beyond peak. */
500 i=peak_i; j=peak_j; local_maxima= -10.0;
501 still_falling=TRUE; lod= 1000.0; starting_value = 2000.0;
502 if(peak_i >= n_intervals) local_minima = 1000.0;
503 else local_minima = interval[i]->point[j]->lod_score;
504 } else { /* Not a local maxima */
505 if (still_falling && lod>prev_lod) still_falling=FALSE;
506 if (!still_falling && lod>local_maxima)
507 { local_maxima=lod; peak_i=i; peak_j=j; }
508 if (lod<local_minima) { local_minima=lod; }
509 }
510 } while (step(&i,&j,&contig,interval,n_intervals,TRUE));
511
512 /* Deal with with fencepost error - have a max if lod is rising at end! */
513 if (!off_end && lod>threshold && lod==local_maxima) {
514 peak=peak_finder(&peak_i,&peak_j,&off_end,get_peak_maps,
515 interval,n_intervals,threshold,confidence_falloff);
516 if (first==NULL) first=last=peak;
517 else { last->next=peak; last=peak; }
518 }
519
520 return(first);
521 }
522
523
peak_finder(interval,point,off_end,get_peak_maps,wig,n_intervals,threshold,falloff)524 WIGGLE_PEAK *peak_finder(interval,point,off_end,get_peak_maps,wig,n_intervals,
525 threshold,falloff)
526 int *interval, *point;
527 bool *off_end;
528 bool get_peak_maps; /* if TRUE, fill in peak->map MAY REQUIRE CALCULATION! */
529 WIGGLE_INTERVAL **wig;
530 int n_intervals;
531 real threshold, falloff;
532 {
533 /* i's are intervals and j's are points */
534 int peak_i, peak_j, left_i, left_j, right_i, right_j, i, j, k;
535 real lod, maxima;
536 bool off_right, off_left, contig;
537 WIGGLE_PEAK *peak=NULL;
538
539 run {
540 *off_end=TRUE;
541
542 /* First go to right boundary, finding the peak in the process. The args
543 will be side-effected to be one step beyond the right boundary. */
544
545 right_i=peak_i= *interval; right_j=peak_j= *point;
546 maxima= wig[*interval]->point[*point]->lod_score;
547 off_right=TRUE; *off_end=TRUE;
548
549 while (step(interval,point,&contig,wig,n_intervals,TRUE)) {
550 if (!contig) { *off_end=FALSE; break; }
551 lod= wig[*interval]->point[*point]->lod_score;
552 if (lod>maxima) { maxima=lod; peak_i= *interval; peak_j= *point; }
553 else if ((falloff<0.0 && lod<maxima+falloff) ||
554 (falloff>0.0 && lod<falloff))
555 { off_right=FALSE; *off_end=FALSE; break; }
556 else { right_i= *interval; right_j= *point; }
557 }
558
559 /* Now go to left boundary... */
560
561 i=left_i=peak_i; j=left_j=peak_j; off_left=TRUE;
562 while (step(&i,&j,&contig,wig,n_intervals,FALSE)) {
563 if (!contig) { break; }
564 lod= wig[i]->point[j]->lod_score;
565 if ((falloff<0.0 && lod<maxima+falloff) ||
566 (falloff>0.0 && lod<falloff)) { off_left=FALSE; break; }
567 else { left_i=i; left_j=j; }
568 }
569
570 /* And finally pack these data away... */
571
572 single(peak,WIGGLE_PEAK);
573 peak->next=NULL; peak->map=NULL;
574 k= wig[peak_i]->map->num_intervals-1;
575
576 peak->left= wig[peak_i]->map->left[k];
577 peak->right= wig[peak_i]->map->right[k];
578 peak->qtl_pos= wig[peak_i]->point[peak_j]->qtl_pos;
579 peak->lod_score= wig[peak_i]->point[peak_j]->lod_score;
580 peak->var_explained= wig[peak_i]->point[peak_j]->var_explained;
581 peak->qtl_weight= wig[peak_i]->point[peak_j]->qtl_weight;
582 peak->qtl_dominance= wig[peak_i]->point[peak_j]->qtl_dominance;
583
584 peak->forward_left= wig[right_i]->map->left[k];
585 peak->forward_right= wig[right_i]->map->right[k];
586 if (off_right) peak->forward_pos=OFF_END;
587 else peak->forward_pos= wig[right_i]->point[right_j]->qtl_pos;
588
589 peak->backward_left= wig[left_i]->map->left[k];
590 peak->backward_right= wig[left_i]->map->right[k];
591 if (off_left) peak->backward_pos=OFF_END;
592 else peak->backward_pos= wig[left_i]->point[left_j]->qtl_pos;
593
594 if (get_peak_maps) {
595 /* We assume here that best_point is pretty much near the
596 ML_qtl_pos, in that we ignore max_like_map. */
597 if (wig[peak_i]->best_point==peak_j) {
598 peak->map=wig[peak_i]->map;
599 peak->free_map=FALSE;
600 } else { /* Shit! no map available for this pos, so calculate it */
601 peak->map= alloc_qtl_map(wig[peak_i]->map->num_intervals,
602 wig[peak_i]->map->num_continuous_vars);
603 peak->free_map=TRUE;
604 mapcpy(peak->map,wig[peak_i]->map);
605 peak->map->fix_pos[k]=peak->qtl_pos;
606 temp_print("Calculating QTL map...");
607 make_qtl_map(peak->map);
608 temp_print(NULL);
609 }
610 }
611
612 } when_aborting { free_wiggle_peaks(peak); relay; }
613 return(peak);
614 }
615
616
free_wiggle_peaks(p)617 void free_wiggle_peaks(p)
618 WIGGLE_PEAK *p;
619 {
620 if (p==NULL) return;
621 free_wiggle_peaks(p->next);
622 if (p->free_map && p->map!=NULL) free_qtl_map(map);
623 unsingle(p,WIGGLE_PEAK);
624 }
625
626
isa_test_wiggle(wiggle_num)627 bool isa_test_wiggle(wiggle_num)
628 int wiggle_num;
629 {
630 QTL_SEQUENCE *p;
631
632 if (raw.data_type==BACKCROSS) return(FALSE);
633 p=wiggles[wiggle_num]->seq; while (p->next!=NULL) p=p->next;
634 if (p->genetics.interx_type!=TEST_MODELS) return(FALSE);
635 return(TRUE);
636 }
637
638
639 #define NO_WIGGLES \
640 "No scan results have yet been saved.\nUse the 'scan' command first."
641 #define WIGGLE_NUM \
642 "%d is not a valid number for a saved scan result.\nUse a number from %d-%d.\n"
643 #define ORDER_ONE \
644 "%d.%d is not a valid number for a saved scan result.\nScan %d only has one entry, numbered %d.1.\n"
645 #define ORDER_NUM \
646 "%d.%d is not a valid number for a saved scan result.\nUse a number from %d.1-%d.%d.\n"
647
get_wiggle_nums(str,wiggle,order)648 void get_wiggle_nums(str,wiggle,order)
649 char *str;
650 int *wiggle, *order;
651 {
652 bool order_set;
653 int i;
654
655 if (num_wiggles==0) error(NO_WIGGLES);
656 if (nullstr(str)) { *wiggle= *order= -1; return; }
657
658 if ((i=strfinder(str,'.'))>0) {
659 str[i++]='\0'; if (sscanf(str,"%d",wiggle)<1) usage_error(num_args);
660 order_set=TRUE; if (sscanf(&str[i],"%d",order)<1)
661 usage_error(num_args);
662 if ((--*wiggle)<0 || (--*order)<0) usage_error(num_args);
663 } else {
664 *order= -1; if (sscanf(str,"%d",wiggle)<1) usage_error(num_args);
665 order_set=FALSE; if ((--*wiggle)<0) usage_error(num_args);
666 }
667
668 if (*wiggle<first_wiggle || *wiggle>=num_wiggles) {
669 sf(ps,WIGGLE_NUM,*wiggle+1,first_wiggle+1,num_wiggles);
670 error(ps);
671 }
672 if (wiggles[*wiggle]==NULL) send(CRASH);
673
674 if (order_set)
675 if (*order>=wiggles[*wiggle]->num_orders) {
676 if (wiggles[*wiggle]->num_orders==1)
677 sf(ps,ORDER_ONE,*wiggle+1,*order+1,*wiggle+1,*wiggle+1);
678 else sf(ps,ORDER_NUM,*wiggle+1,*order+1,*wiggle+1,*wiggle+1,
679 wiggles[*wiggle]->num_orders);
680 error(ps);
681 }
682 }
683
684
685 #define NO_COMPARES \
686 "No compare results have yet been saved.\nUse the 'compare' command first."
687 #define COMPARE_NUM \
688 "%d is not a valid number for a saved compare result.\nUse a number from %d-%d.\n"
689 #define ORDER_ONE_COMP \
690 "%d.%d is not a valid number for a saved compare result.\nCompare %d only has one contiguous entry, numbered %d.1.\n"
691 #define ORDER_NUM_COMP \
692 "%d.%d is not a valid number for a saved compare result.\nUse a number from %d.1-%d.%d.\n"
693
694
get_compare_nums(str,compare,contig)695 void get_compare_nums(str,compare,contig)
696 char *str;
697 int *compare,*contig;
698 {
699 int i;
700 bool contig_set;
701
702 if(num_compares == 0) error(NO_COMPARES);
703 if (nullstr(str)) { *compare= *contig= -1; return; }
704
705 if((i=strfinder(str,'.')) > 0) {
706 str[i++] = '\0'; if(sscanf(str,"%d",compare) < 1) usage_error(num_args);
707 contig_set = TRUE;
708 if(sscanf(&str[i],"%d",contig)<1) usage_error(num_args);
709 if ((--*compare)<0 || (--*contig)<0) usage_error(num_args);
710 } else {
711 *contig = -1; if(sscanf(str,"%d",compare)<1) usage_error(num_args);
712 contig_set = FALSE; if ((--*compare)<0) usage_error(num_args);
713 }
714
715 if(*compare == first_compare || *compare >= num_compares) {
716 sf(ps,COMPARE_NUM,*compare+1,first_compare+1,num_compares);
717 error(ps);
718 }
719 if(compares[*compare]==NULL) send(CRASH);
720
721 if(contig_set)
722 if(*contig>=compares[*compare]->num_contigs) {
723 if(compares[*compare]->num_contigs==1)
724 sf(ps,ORDER_ONE_COMP,*compare+1,*contig+1,*compare+1,*compare+1);
725 else
726 sf(ps,ORDER_NUM_COMP,*compare+1,*contig+1,*compare+1,*compare+1,
727 compares[*compare]->num_contigs);
728 error(ps);
729 }
730 }
731
732
save_wiggle(fp,n)733 void save_wiggle(fp,n)
734 FILE *fp;
735 int n;
736 {
737 int i,j;
738 if (wiggles[n]->data == NULL)
739 fprintf(fp,"%d\n",-1);
740 else {
741 fprintf(fp,"%d %s\n",wiggles[n]->trait,wiggles[n]->seq_string);
742
743 for(i = 0; i < wiggles[n]->num_orders; i++) {
744 for(j = 0; j < wiggles[n]->num_wiggled_intervals; j++) {
745 if(wiggles[n]->data[i][j] != NULL) {
746 save_interval(fp, wiggles[n]->data[i][j]);
747 }
748 }
749 }
750 }
751 }
752
753
save_interval(fp,interval)754 void save_interval(fp,interval)
755 FILE *fp;
756 WIGGLE_INTERVAL *interval;
757 {
758 int i;
759
760
761 sprintf(ps,"%d %.3lf %d %d %d %d\n",interval->num_points,
762 interval->cm_increment, interval->contig, interval->max_like_map,
763 interval->in_qtls_struct, interval->best_point);
764 fpr(fp);
765
766 /* save the peak map */
767 save_qtl_map(fp,interval->map);
768
769 /* save other data points (if any) */
770 for(i = 0; i < interval->num_points; i++) {
771 sf(ps,"%.4lf %.3lf %.3lf %.3lf %.4lf\n",interval->point[i]->qtl_pos,
772 interval->point[i]->qtl_weight, interval->point[i]->qtl_dominance,
773 interval->point[i]->lod_score, interval->point[i]->var_explained);
774 fpr(fp);
775 }
776 }
777
778
save_qtl_map(fp,map)779 void save_qtl_map(fp,map)
780 FILE *fp;
781 QTL_MAP *map;
782 {
783 int i;
784
785 sf(ps,"%d %d %d %d %d\n",map->trait,map->num_intervals,map->max_intervals,
786 map->num_continuous_vars,map->max_continuous_vars);
787 fpr(fp);
788
789 for(i = 0; i < map->num_intervals; i++) {
790 sf(ps,"%d %d %.3lf %.3lf %.3lf %.3lf %.3lf\n", map->left[i],
791 map->right[i], map->interval_len[i], map->qtl_pos[i],
792 map->qtl_weight[i], map->qtl_dominance[i], map->fix_pos[i]);
793 fpr(fp);
794 sf(ps,"%.4lf %d %.4lf %.4lf %.4lf\n", map->constraint[i].backx_weight,
795 map->constraint[i].interx_type, map->constraint[i].a,
796 map->constraint[i].b, map->constraint[i].c);
797 fpr(fp);
798 }
799
800 for(i = 0; i < map->num_continuous_vars; i++) {
801 sf(ps,"%d %lf %lf\n",map->cont_var[i],map->cont_var_weight[i],
802 map->fix_cont_var_weight[i]);
803 fpr(fp);
804 }
805
806
807 sf(ps,"%.4lf %.4lf %.4lf %.4lf %.4lf %.4lf %.4lf %.4lf %.4lf %.4lf\n",
808 map->mu, map->sigma_sq, map->var_explained, map->chi_sq,
809 map->null_mu, map->null_sigma_sq, map->null_log_like,
810 map->log_like, map->no_data_like, map->abs_log_like);
811 fpr(fp);
812 }
813
814
load_wiggle(fp)815 void load_wiggle(fp)
816 FILE *fp;
817 {
818 int trait,i,j,n;
819 int n_orders,n_ints,n_wiggled,foo;
820 QTL_SEQUENCE *seqnce;
821
822 fgetdataln(fp,NULL);
823 if(!itoken(&ln,iREQUIRED,&trait)) return;
824 if (trait != -1) {
825 _filter(ln);
826 seqnce = compile_intervals(ln);
827 reset_state(seqnce,TRUE,&n_ints,&foo,&n_orders,&n_wiggled);
828 }
829 n = allocate_wiggle_struct(trait,seqnce,ln,n_ints,n_orders,n_wiggled);
830
831 if(n == -1) return;
832 if (trait != -1) {
833 for(i = 0; i < wiggles[n]->num_orders; i++) {
834 for(j = 0; j < wiggles[n]->num_wiggled_intervals; j++) {
835 if(wiggles[n]->data[i][j] != NULL) {
836 load_interval(fp, wiggles[n]->data[i][j]);
837 }
838 }
839 }
840 }
841 }
842
843
load_interval(fp,interval)844 void load_interval(fp,interval)
845 FILE *fp;
846 WIGGLE_INTERVAL *interval;
847 {
848 int num_points,i,matched,i1,i2,i3,i4;
849 real pos,weight,dom,lod,var_exp,inc;
850
851 run {
852 if(fscanf(fp,"%d %lf %d %d %d %d\n",&num_points,&inc,&i1,&i2,&i3,&i4)
853 != 6)
854 error("error in loading saved qtls file");
855
856 interval->num_points = num_points; interval->cm_increment = inc;
857 interval->contig = i1; interval->max_like_map = i2;
858 interval->in_qtls_struct = i3; interval->best_point = i4;
859
860 /* load map */
861 interval->map = load_qtl_map(fp);
862
863 /* load other data points */
864 if(interval->num_points > 0)
865 parray(interval->point, interval->num_points, WIGGLE_POINT);
866
867 for(i = 0; i < interval->num_points; i++) {
868 matched = fscanf(fp,"%lf %lf %lf %lf %lf\n",
869 &pos,&weight,&dom,&lod,&var_exp);
870
871 if(matched != 5)
872 error("error in loading saved qtls file");
873
874 interval->point[i]->qtl_pos = pos;
875 interval->point[i]->qtl_weight = weight;
876 interval->point[i]->qtl_dominance = dom;
877 interval->point[i]->lod_score = lod;
878 interval->point[i]->var_explained = var_exp;
879 }
880 } when_aborting {
881 if(msg == ENDOFILE) {
882 print("Error while loading saved qtls file...data not loaded.\n");
883 }
884 else relay;
885 }
886 }
887
888
load_qtl_map(fp)889 QTL_MAP *load_qtl_map(fp)
890 FILE *fp;
891 {
892 int i1,i2,i3,i4,i5, i;
893 real r1,r2,r3,r4,r5,r6,r7,r8,r9,r10;
894 QTL_MAP *map;
895
896 if(fscanf(fp,"%d %d %d %d %d\n",&i1,&i2,&i3,&i4,&i5) != 5)
897 error("error in loading qtls file");
898
899 map = alloc_qtl_map(i3,i5);
900
901 map->trait = i1; map->num_intervals = i2;
902 map->max_intervals = i3;
903 map->num_continuous_vars = i4; map->max_continuous_vars = i5;
904
905 for(i = 0; i < map->num_intervals; i++) {
906 if(fscanf(fp,"%d %d %lf %lf %lf %lf %lf\n",&i1,&i2,&r1,&r2,&r3,&r4,&r5)
907 != 7)
908 error("error in loading qtls file");
909
910 map->left[i] = i1; map->right[i] = i2;
911 map->interval_len[i] = r1; map->qtl_pos[i] = r2;
912 map->qtl_weight[i] = r3; map->qtl_dominance[i] = r4;
913 map->fix_pos[i] = r5;
914
915 if(fscanf(fp,"%lf %d %lf %lf %lf\n",&r1,&i1,&r2,&r3,&r4) != 5)
916 error("error in loading qtls file");
917
918 map->constraint[i].backx_weight = r1;
919 map->constraint[i].interx_type = i1; map->constraint[i].a = r2;
920 map->constraint[i].b = r3; map->constraint[i].c = r4;
921 }
922
923 for(i = 0; i < map->num_continuous_vars; i++) {
924 if(fscanf(fp,"%d %lf %lf\n",&i1,&r1,&r2) != 3)
925 error("error loading qtls file");
926
927 map->cont_var[i] = i1; map->cont_var_weight[i] = r1;
928 map->fix_cont_var_weight[i] = r2;
929 }
930
931 if(fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
932 &r1,&r2,&r3,&r4,&r5,&r6,&r7,&r8,&r9,&r10) != 10)
933 error("error in loading qtls file");
934
935 map->mu = r1; map->sigma_sq = r2;
936 map->var_explained = r3; map->chi_sq = r4;
937 map->null_mu = r5; map->null_sigma_sq = r6;
938 map->null_log_like = r7; map->log_like = r8;
939 map->no_data_like = r9; map->abs_log_like = r10;
940
941 return(map);
942 }
943
944
save_compare(fp,n)945 void save_compare(fp,n)
946 FILE *fp;
947 int n;
948 {
949 int i;
950
951 if(compares[n] == NULL) return;
952 if (compares[n]->data == NULL) {
953 fprintf(fp,"%d\n", -1);
954 return;
955 }
956 fprintf(fp,"%d %s\n",compares[n]->trait,compares[n]->seq_string);
957
958 for(i = 0; i < compares[n]->num_orders; i++) {
959 if(compares[n]->data[i] != NULL) {
960 fprintf(fp,"%d\n",compares[n]->data[i]->contig);
961 save_qtl_map(fp,compares[n]->data[i]->map);
962 }
963 }
964 }
965
966
load_compare(fp)967 void load_compare(fp)
968 FILE *fp;
969 {
970 int trait,i,contig,n;
971 int n_orders,n_ints,n_wiggled,foo;
972 QTL_SEQUENCE *seqnce;
973
974 fgetdataln(fp,NULL);
975 if(!itoken(&ln,iREQUIRED,&trait)) return;
976 if(trait != -1) {
977 _filter(ln);
978 seqnce = compile_intervals(ln);
979 reset_state(seqnce,FALSE,&n_ints,&foo,&n_orders,&n_wiggled);
980 }
981 n = allocate_compare_struct(trait,seqnce,ln,n_ints,n_orders);
982 if(n == -1) return;
983 if(compares[n]->num_contigs != -1) {
984 compares[n]->num_contigs = 1;
985 for(i = 0; i < compares[n]->num_orders; i++) {
986 fscanf(fp,"%d\n",&contig);
987 compares[n]->data[i]->contig = contig;
988 if(!contig) compares[n]->num_contigs += 1;
989 compares[n]->data[i]->map = load_qtl_map(fp);
990 }
991 }
992 }
993
994
995
996
997