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 /***** MAP struct - holds all information about a linkage map - maps.c *****/ 15 typedef struct { 16 char *map_name; /* name of MAP */ 17 int max_loci; 18 int num_loci; 19 int *locus; /* [num_loci] (markers in the map) */ 20 double **rec_frac; /* [intervals][M/F], recombination fractions */ 21 int unlink; /* indicates interval (if any) held unlinked */ 22 int *fix_interval; /* [intervals], intervals not converging */ 23 int sex_specific; 24 double log_like; 25 bool allow_errors; /* if TRUE, below things are valid */ 26 double **error_lod; /* [num_loci][indiv] - maybe NULL */ 27 double *error_rate; /* [num_loci] - NULL if error_lod==NULL */ 28 bool modified; 29 } MAP; 30 31 /* for map->unlink */ 32 #define MANY_UNLINKED (-2) 33 #define NONE_UNLINKED (-1) 34 /* A number 1 or greater indicates the interval to be held unlinked (interval 35 0 is the open interval left of the right-most locus). */ 36 37 /* Functions for MAP struct in maps.c */ 38 MAP *allocate_map(); /* arg: max_loci; alloc map for up to max_loci */ 39 void free_map(); /* arg: MAP *map; frees map previously allocated */ 40 void mapcpy(); /* args: MAP *to, *from; bool clean_it; copy from->to */ 41 int insert_locus(); /* args: map, position, locus; adds locus to map */ 42 bool clean_map(); /* arg: map; resets rf's etc in map, returns TRUE */ 43 void setup_error_matrix(); /* args: MAP *map, int n_indivs, real rate; */ 44 45 /* Call init_rec_fracs only if you know that map->rec_frac[][] have been 46 initialized to either NOT_FIXED, UNLINK_ME, or a fixed value. */ 47 void init_rec_fracs(); /* arg: MAP *m; sets rf, fix_interval, etc for ctm */ 48 void init_not_fixed(); /* arg: MAP *m; sets rf, fix_interval for ctm */ 49 void init_for_ctm(); /* arg: MAP *m; bool sex, errors, hot_start; similar */ 50 51 52 /***** SAVED_LIST Multiple map structure - maps.c *****/ 53 54 typedef struct { 55 int max_maps; /* used in allocation of structure */ 56 int max_loci; /* used in allocation of structure */ 57 MAP **map_list; /* list of max_maps pointers to MAPs */ 58 MAP *extra_map; /* expendable map-for insertion and deletion */ 59 int num_maps; /* number of maps in map_list (<= max_maps) */ 60 double threshold; /* highest acceptable likelihood when sorting */ 61 bool sorted; /* determines whether this is a sorted list */ 62 } SAVED_LIST; 63 64 #define VERY_UNLIKELY -999999.0 /* useful threshold */ 65 #define SORTED 1 66 #define UNSORTED 0 67 #define FULL_LIST (-1) 68 69 SAVED_LIST *allocate_map_list(); /* args: max_maps, max_loci, sorted; */ 70 MAP *get_map_to_bash(); /* arg: SAVED_LIST *list; gets the extra map to bash */ 71 int insert_map_into_list(); /* arg: *SAVED_LIST, **MAP; inserts *MAP */ 72 /* sets *MAP to new extra_map for re-use */ 73 void overwrite_map_num(); /* args: SAVED_LIST *list; MAP **map; int chrom; */ 74 75 void free_map_list(); /* arg: *SAVED_LIST; frees SAVED_LIST structure */ 76 void clean_list(); /* arg: *SAVED_LIST; resets values in list for re-use */ 77 MAP *get_best_map(); /* arg: *SAVED_LIST; returns best map in list */ 78 79 void read_map(), write_map(); /* args: FILE *fp; */ 80 81 82 83 /***** Global structs used to store 2pt and 3pt data - in map_23.c *****/ 84 85 typedef struct { 86 double lodscore[3]; /* actually, only SEXSPEC and NOSEX are used */ 87 double theta[3]; /* MALE, FEMALE, and NOSEX */ 88 bool used; /* for garbage collection */ 89 } TWO_PT_DATA; 90 91 #define UNLINKED_LOD 0.50 /* Because LODs less than this never get stored */ 92 #define UNLINKED_THETA 0.40 /* we should not allow lodbounds<UNLINKED_LOD */ 93 94 void allocate_two_pt(); /* arg: num_loci; allocs and inits global 2pt data */ 95 void free_two_pt(); /* arg: num_loci; frees global */ 96 void expand_two_pt(); /* arg: num_entries; pre-allocates space in global */ 97 #define EXPAND_DEFAULT (-1) /* can be arg to above */ 98 99 /* these two calculate on demand WHAT ABOUT HAPS??? */ 100 void compute_two_pt(); /* args: int locus1, locus2; */ 101 bool get_two_pt(); /* args: int locus1, locus2; real *lod, *theta; */ 102 /* Returns TRUE if !unlinked */ 103 bool get_sex_two_pt(); /* args: int locus1, locus2; real *lod, *m, *f; */ 104 105 void allocate_three_pt(); /* arg: num_loci; allocs and inits global 3pt data */ 106 void free_three_pt(); /* arg: num_loci; frees global */ 107 void bash_all_three_pt(); /* arg: num_loci; frees then allocates */ 108 109 bool restore_triple(); 110 /* args: int locus1, locus2, locus3; real *delta1, *delta2, *delta3; 111 returns TRUE and side-effects deltas if three-pt data exist */ 112 113 void compute_triple(); 114 /* args: int locus1, locus2, locus3; real *delta1, *delta2, *delta3; 115 log-likelihoods (0.0 or negative) are filled in. */ 116 117 void insert_triple(); 118 /* args: int locus1, locus2, locus3; real delta1, delta2, delta3; saves it */ 119 120 bool three_linked(); /* args: int *locus; real lod, theta; int nlinks; */ 121 void compute_3pt(); /* args: *seq; bool sex; real error_rate, *like; *map; */ 122 123 extern bool two_pt_touched; 124 extern bool three_pt_touched; 125 126 void read_two_pt(); 127 void write_two_pt(); 128 void read_three_pt(); 129 void write_three_pt(); 130 131 132 /**************** Order building code in orders.c ****************/ 133 134 void get_linkage_group(); 135 136 void setup_haplo_group(); /* args: loci, num_loci; */ 137 bool delete_haplo_groups(); /* args: loci, num_loci; TRUE if any deleted */ 138 void find_haplo_group(); /* args: loci *#loci group *#group obs1 obs2 */ 139 bool force_haplo_sanity(); /* int *locus; bool verbose; FALSE if changed */ 140 141 void setup_3pt_data(); 142 /* args: int locus[], num_loci; real three_pt_threshold; 143 Sets the global matrix used by create_order. three_pt_threshold 144 should be negative, or if zero, 3pt data are unused. */ 145 void free_3pt_data(); 146 147 bool start_3pt(); 148 /* args: int locus[], *num_loci; real threshold; int order[], *num_ordered; 149 If TRUE returned, all args (but threshold) are side-effected. */ 150 151 int three_pt_exclusions(); 152 /* args: int order[], num_placed, new_locus; bool excluded[]; 153 returns num_orders; 154 Examines three-pt data for inserting new_locus in a order. 155 excluded[i], for i=0...num_placed, and *num_orders are side-effected. 156 This assumes that excluded[i] is initialized (likely to FALSE) */ 157 158 bool three_pt_verify(); 159 /* args: int loci[], num_loci, window_size; 160 Examines an order to see if it is compatible with the three-point data. 161 Return FALSE if order is excluded. */ 162 163 void informative_subset(); 164 /* args: locus[], num_loci, min_infs; real min_dist; bool codom_only, haplos; 165 int subset[], *num; finds a informative and not-too-closely-spaced set 166 of loci. Does not understand sex_specific, nor does it really do anything 167 smart with mixed typings. */ 168 169 void extend_order(); 170 /* int create_order(); OBSOLETE */ 171 /* args: int *loci, *num_loci, *order, *num_ordered; real npt_threshold; 172 Iteratively tries to place loci in order - side-effecting order[i] and 173 *num_ordered. Is rather chatty. */ 174 175 void find_window(); 176 /* int order[], num_placed, new_locus, excluded[], min_window, *start, *end; 177 find the appropriate region in which to try to place a marker, based on 178 excluded[], and maybe (someday) informativeness */ 179 180 void place_locus(); 181 /* args: MAP *map; int locus; int start, finish; bool excluded[]; 182 PLACE *result[]. int *best_pos; MAP *best_map, *temp_map; 183 Tries locus into map order between locus #s i=start..finish 184 (inclusive), for which !excluded[i]. result[i] is side-effected, where 185 like may be 0.0 (best), NO_LIKE (untried) or ZERO_LIKE (if places on 186 top of a marker, one side gets this value */ 187 188 int npt_exclusions(); /* return num_ok_intervals */ 189 /* args: PLACE *placements[]; int num_loci; real like_threshold; 190 real worst_error_threshold, net_error_threshold; 191 bool excluded[], *zero_placement, *placed_off_end; real *error_lod; 192 bool *single_error; real *second_best_like, *best_unallowed_like; 193 Examines the likes produced by place_locus(), side-effecting excluded[i], 194 and the many flags. */ 195 196 typedef struct { /* nifty little struct used for the above */ 197 real like; 198 real dist; 199 bool zero; 200 real net_error; 201 real worst_error; 202 } PLACE; 203 204 typedef struct { /* working data for place cmd and extend_order() */ 205 int locus; 206 bool *excluded; /* [position-in-order] */ 207 int num_places, best_pos; 208 bool off_end, errors; 209 int priority; /* a random number */ 210 MAP *best_map; 211 } PLACEME; 212 213 214 215 /***** Chromosome framework, assignment, and placement stuff - chroms.c *****/ 216 217 extern SAVED_LIST *chromosome; /* malloced by allocate_mapping_structs() */ 218 #define chrom2str(x) ((x)>=0 ? chromosome->map_list[x]->map_name : WRS("none")) 219 bool isa_chrom(); /* args: char *name; int *chrom; side-effected if TRUE */ 220 #define num_chromosomes (chromosome->num_maps) 221 extern int current_chrom; /* set by the sequence command or reset_state() */ 222 223 bool make_new_chrom(); /* args: char *name; returns # */ 224 void set_chrom_frame(); /* int chrom; char *name; MAP *map; */ 225 /* When this is called, new MUST be equal to get_map_to_bash(chromosome), amd 226 assigned_to(*,chrom) must have be true for all loci in new map. We also 227 assume haplo_sanity is true for old framework, and force it for new one. */ 228 229 MAP *get_chrom_frame(); /* args: int chrom, *n_loci; n_loci may be 0 or 1 */ 230 bool framework_marker(); /* int locus; haplos are NOT checked: if use_haplos 231 is on, the chrom frames better only include the first marker in any haplo 232 group */ 233 234 void get_chrom_loci(); /* args: int chrom, *locus, *num_loci; int which_loci; 235 copies the chrom's loci into an array. */ 236 #define FRAMEWORK 1 /* do not change these w/o changing get_chrom_loci() */ 237 #define NON_FRAME 2 238 #define ALL_LOCI 3 /* =FRAMEWORK+NON_FRAME */ 239 #define HAPLOS 4 240 /* Rule on haplos: if use_haplotypes is on then only the first locus in any 241 haplo group is included, unless HAPLOS is specified. if use_haplotypes is, 242 off then all fw/assigned/placed markers are included. */ 243 244 void count_chrom_loci(); 245 /* args: int chrom, *n_anchor, *n_frame, *n_total, *n_placed, *n_unique; 246 int *temp; bool haplos_too; temp should be aloced raw.markers long */ 247 248 typedef struct { 249 int chromosome; /* framework this marker is linked to */ 250 int linked_to; /* closest linked marker in that framework */ 251 real LODscore, theta; /* two-point values with that marker */ 252 bool modified; /* needs to be saved */ 253 int status; 254 } ASSIGNMENT; 255 256 extern ASSIGNMENT **assignment; 257 258 /* assignment status values: */ 259 /* Data which must be valid: CHROM LOD&LOCUS NAME */ 260 #define A_PROBLEM (-3) /* nyet nyet conflict */ 261 #define A_CHANGED (-2) /* nyet nyet changed */ 262 #define A_UNLINKED (-1) /* nyet nyet unlinked */ 263 #define A_UNKNOWN (0) /* nyet nyet - */ 264 #define A_BORDERLINE (1) /* da da low-lod */ 265 #define A_ASSIGNED (2) /* da da linked */ 266 #define A_ATTACH (3) /* da nyet attached */ 267 #define A_FRAMEWORK (4) /* da nyet framework */ 268 #define A_ANCHOR (5) /* da nyet anchor */ 269 /* note: A_FRAMEWORK means that the locus was in a framework when 270 do_assignment was called, resulting in loss of linkage info for it */ 271 272 /* all four below silently force_haplo_sanity */ 273 bool assigned(); /* args: int locus; */ 274 bool assigned_to(); /* args: int locus, chrom; TRUE if locus is on chrom */ 275 bool anchor_locus(); /* int locus; */ 276 int assignment_chrom(); /* args: int locus; */ 277 int assignment_state(); /* args: int locus; */ 278 279 /* Macros for yucks, implicitly require force_haplo_sanity */ 280 #define assignment_lod(locus) (assignment[locus]->LODscore) 281 #define assignment_locus(locus) (assignment[locus]->linked_to) 282 #define assignment_theta(locus) (assignment[locus]->theta) 283 284 bool is_assignable(); /* args: int locus, chrom; prints msg if FALSE */ 285 /* note that is_assignable is FALSE if locus is insane, haplo-wise */ 286 287 void assign_this(); 288 /* int locus, state, chrom; real lod, theta; int linked_locus; char *msg; 289 msg pre-empts any other message, only used for A_PROBLEM as yet 290 This can unassign anchor markers, but not reassign them! */ 291 292 void attach_this(); /* int locus, state; char *msg; just a convenience */ 293 void unassign_this(); /* int locus, state; char *msg; just a convenience */ 294 295 void set_chrom_anchors(); /* int chrom, *locus, num_loci */ 296 /* The only right way to assign() a locus as A_ANCHOR. Will reassign if need 297 be and it can. Do not use is_assignable() beforehand, as that never lets 298 one assign an anchor. */ 299 300 void do_assignments(); /* args: int *locus, num_loci; 301 real lod1, unlinked_lod1, theta1, lod2, unlinked_lod2, theta2; bool haplo; 302 Does the real work - is_assignable must have been verified */ 303 304 typedef struct { 305 int chromosome; 306 int num_intervals; 307 int interval[MAX_STORED_PLACES]; 308 real like_ratio[MAX_STORED_PLACES]; 309 real distance[MAX_STORED_PLACES]; 310 real errors[MAX_STORED_PLACES]; 311 real threshold; 312 real error_lod; 313 bool single_error; 314 bool modified; 315 int status; 316 } PLACEMENT; 317 318 extern PLACEMENT **placement; 319 320 /* placement status values: */ 321 /* Data which must be valid: LIKES ERRORS NAME */ 322 #define M_PROBLEM (-2) /* nyet nyet problem */ 323 #define M_FRAMEWORK (-1) /* nyet nyet framework */ 324 #define M_UNKNOWN (0) /* nyet nyet - */ 325 #define M_ERROR (1) /* da da errors? */ 326 #define M_OFFEND (2) /* da da off-end */ 327 #define M_REGION (3) /* da da region */ 328 #define M_UNIQUE (4) /* da da unique */ 329 #define M_ZERO (5) /* da da zero */ 330 331 #define NO_CHROM (-1) 332 #define ANY_CHROM (-2) /* flag for is_assignable() */ 333 #define NEW_CHROM (-3) /* flag for set_chrom_frame */ 334 #define NO_LOCUS (-1) 335 #define NO_LOD 0.0 336 #define NO_THETA 0.4996 /* HMM_MAX_THETA 0.4995 */ 337 #define NO_LIKE (-999.98765) /* hope that's obscure */ 338 #define ZERO_LIKE (-999.00123) /* hope that's obscure too */ 339 #define NO_ERRORS NO_LIKE /* for ->error_lod */ 340 341 #define FRAMEWORK_MARKER (-1) 342 343 bool placed_locus(); /* args: int locus; forces haplo_sanity w/no msg */ 344 int placement_state(); /* args: int locus; forces haplo_sanity w/no msg */ 345 346 bool is_placeable(); /* args: int locus, new_chrom; */ 347 /* prints msg if FALSE, and it WILL be FALSE if !haplo_sanity */ 348 349 int place_this(); 350 /* args: int locus, chrom; PLACE **placements; real like_threshold; 351 real worst_error_threshold, net_error_threshold; bool excluded[]; 352 sets placement struct and returns #intervals within threshold. 353 uses npt_exclusions() to do its real work. Is chatty. assumes 354 haplo_sanity */ 355 void unplace_this(); /* args: int locus, chrom; bool problem, print_msg; */ 356 357 int best_placement(); /* args: int locus; */ 358 int second_best_placement(); /* args: int locus; real *like; */ 359 360 void allocate_mapping_data(), free_mapping_data(); /* args: num_markers; */ 361 void write_mapping_data(), read_mapping_data(); /* args: FILE *fp; */ 362 void bash_mapping_data(); /*args: int changed_markers[], num_changed; */ 363 364 void print_ps_map(); /* args: fp, map; in ps_maps.c */ 365 void print_ps_chrom(); /* args: fp, chrom; in ps_maps.c */ 366 void print_all_ps_chroms(); /* args: fp; in ps_maps.c */ 367 368 369 /***************** Misc Saved Mapping Info *****************/ 370 /* support is all now in map_123.c */ 371 372 extern char **note; /* [locus]=>string, Note about a particular marker */ 373 extern bool *modified; /* [locus] FALSE if no, TRUE if yes */ 374 extern real *error_rate; /* [locus] apriori, see use_error_rate, state.c */ 375 extern int *my_group; /* [locus] For group command, NO_GROUP if unknown */ 376 extern int num_groups; 377 #define NO_GROUP (-1) 378 379 extern int *haplo_first, *haplo_next; /* [locus] */ 380 extern int *order_first, *unorder_first; /* [order]=> locus# or NO_LOCUS */ 381 extern int *order_next; /* [locus]=>locus or NO_LOCUS, for order or unorder */ 382 extern int num_orders; 383 #define haplotyped(locus) (haplo_first[locus]!=NO_LOCUS) 384 #define haplotype_subordinate(n) (haplo_first[n]!=NO_LOCUS &&haplo_first[n]!=n) 385 386 extern int *class; /* [locus] class of each marker */ 387 extern char **class_name; /* [NUM_CLASSES][NAME_LENGTH] name of each class */ 388 #define NO_CLASS 0 389 #define NUM_CLASSES 11 390 391 bool isa_class(); /* args: char *name; int *num; */ 392 bool make_new_class(); /* args: char *name; char **why_not; */ 393 void print_class_names(); /* lets print() auto_wrap, no nl on end */ 394 395 void allocate_order_data(); /* args: num_markers; */ 396 void free_order_data(); /* args: num_markers; */ 397 void write_order_data(); /* args: FILE *fp; */ 398 void read_order_data(); /* args: FILE *fp; */ 399 void bash_order_info(); /* args: int changed_marker[], num_changed; */ 400 401 402 /*************** Haldane/Kosambi Map Functions - now in maps.c ***************/ 403 typedef struct { 404 char name[12]; 405 double (*add)(); 406 double (*apportion)(); 407 double (*rec_to_dist)(); 408 double (*dist_to_rec)(); 409 double (*d_to_r_deriv)(); 410 } MAP_FUNCTION; 411 412 #define TEN_CM 0.0906346 /* haldane, that is */ 413 #define THIRTY_CM 0.2255942 /* ditto */ 414 #define FIFTY_CM 0.3160603 /* ditto */ 415 #define EIGHTY_CM 0.3990517 /* ditto */ 416 #define ONE_CM 0.00900663 417 #define UNLINKED_CM 0.4995441 /* ditto 350 cM */ 418 #define cm_dist(theta) ((*mapfunction->rec_to_dist)(theta)*100.0) 419 420 extern MAP_FUNCTION *mapfunction; 421 extern MAP_FUNCTION maps[2]; 422 extern int num_map_functions; 423 void map_func(); /* args: int num; sets map function */ 424 #define map_func_name() (mapfunction->name) 425 #define map_func_num() (mapfunction->name[0]=='H' ? HALDANE:KOSAMBI) 426 #define HALDANE 0 427 #define KOSAMBI 1 428