Documentation for CRI-MAP, version 2.4 (3/26/90)
Phil Green, Kathy Falls, and Steve Crooks
INTRODUCTION The main purpose of CRI-MAP is to allow rapid, largely automated construction of multilocus linkage maps (and facilitate the attendant tasks of assessing support relative to alternative locus orders, generating LOD tables, and detecting data errors). Although originally designed to handle codominant loci (e.g. RFLPs) scored on pedigrees "without missing individuals", such as CEPH or nuclear families, it can now (with some caveats described below) be used on general pedigrees, and some disease loci. This version of CRI-MAP is distributed as source code in the language C, on a 360K DOS formatted diskette. To use it, you will need to obtain access to a C compiler for your computer, and compile the program yourself (see the section GETTING STARTED, below). The present version of the program adheres fairly rigidly to the conventions of "standard C" as described in Kernighan and Ritchie (1978) (the main exception being that some variable and function names have more than 8 significant characters), and should be compatible with most C compilers. The program was developed and tested using vers. 2.0 and 2.3 of the VAX C compiler, on a Microvax II with 5 Mb memory under the MicroVMS operating system. It has been successfully ported to a number of other computers, including other VAXes, SUN and Apollo workstations, and the Mac II. CRI-MAP requires a lot of memory; it is desirable to run it on a computer with at least 1 Mb (RAM or virtual memory) if you will be analyzing more than 10 loci simultaneously. It may be possible to run it on an IBM AT under DOS if your data set is small and you reduce the default memory allocations, although we have not tried this yet. A small data set with several chromosome 7 markers is provided with the program for the purpose of testing it (only). I would appreciate being informed of any difficulties in implementing the program, bugs, errors or gaps in the documentation, or suggestions for improvement. Historical note: Version 1 of CRI-MAP was originally written by me in the summer of 1986 in the language APL; the portion of that version which does maximum likelihood estimation for a fixed locus order was based on algorithms developed in collaboration with Eric Lander (Lander and Green, 1987). Collaborative's chromosome 7 map (Barker et al., 1987) was constructed using that version of CRI-MAP, running on an IBM XT. In the summer of 1987 parts of the original version were translated into C, with the help of Steve Crooks, and used in constructing the genome map published in Cell (Eric and his group at MIT independently constructed maps using the program MAPMAKER). At this time I also discovered the "layered EM" maximum likelihood search method, described in (Green, 1988). In January, 1988, I worked out a much faster algorithm for likelihood calculation (the method of switch algebras -- actually a family of related algorithms), also described in (Green, 1988). I have since written and tested the code for these new algorithms and incorporated them into version 2.0 of the program; Kathy Falls has worked on improving various other parts of the program. Versions subsequent to 2.3 were developed by me at Washington University. I soon hope to incorporate a full likelihood analysis for pedigrees with missing individuals, and allow for disease loci with incomplete penetrance. Any user feedback will be gratefully received, and I will do my best to incorporate suggestions into future versions of the program. Phil Green Dept. of Genetics Box 8232 Washington University School of Medicine St. Louis, MO 63110 (314) 362-5192 FAX: (314) 362-4137 EMAIL: pg@genome.wustl.edu CONTENTS 1. Using CRI-MAP with general pedigrees and disease loci 2. Getting started 3. File structures .gen file .par file .dat file .ord file 4. Program options all build chrompic fixed flipsn instant merge prepare quick twopoint 5. Technical notes 6. Changes incorporated in versions 2.2 - 2.4 7. Bibliography USING CRI-MAP WITH GENERAL PEDIGREES AND DISEASE LOCI CRI-MAP can only handle disease loci to the extent that genotypes are known: this means that with autosomal dominant loci, affected individuals are to be scored as heterozygotes at the disease locus, while unaffecteds are scored either as homozygotes for the normal allele (if you are willing to assume complete penetrance for that individual), or as genotype missing (alleles both 0) otherwise. For recessive loci, only affecteds and obligate carriers should be scored for the disease locus. With general pedigrees, CRI-MAP deduces missing genotypes where possible, and computes a likelihood based only on analysis of the known or deduced genotypes. (Partial genotypes, i.e. those in which only one allele is known or deducible, are also utilized in the likelihood calculation to the extent possible). This likelihood is the usual likelihood (as defined for example by Ott (1985)) when the pedigree data are "complete" (i.e. genotypes are known or can be deduced at every locus, for every individual having descendants in the pedigree). If data are missing for an individual at a particular locus, and the possible genotypes at that locus (i.e. those genotypes compatible with the genotypes of ancestors and descendants) include a homozygous genotype, then all meioses in that individual are treated by CRI-MAP as uninformative for that locus. A full likelihood analysis (as provided for example by LIPED or LINKAGE) would instead consider in turn each possible genotype at the locus, assign relative probabilities to each based on the genotypes of ancestors and descendants, and compute a likelihood which is the weighted sum of the likelihood expressions for each particular choice of genotype. CRI-MAP thus ignores some of the available information. In cases that we have examined, however, the information loss appears to be small. If the missing locus genotype is in an "original parent" (i.e. an individual with no ancestors in the pedigree), then in a full likelihood analysis the population allele frequencies are used to assign probabilities to various possible genotypes. These enter into the likelihood by influencing the probability that the allele in any child of the original parent is derived from that parent. CRI-MAP does not make use of allele frequencies; for any allele in a child of an original parent, CRI-MAP determines the parental origin when this can be deduced from the other genotype information, but otherwise assigns equal probability to the two possible parental origins for the allele. In fact, little information is lost by this procedure, except when the allele is rare. For a rare dominant disease, if one parent is known to be affected, the other should therefore be scored as unaffected (rather than as missing, for example). It should be noted in any case that the frequencies of RFLP alleles have usually been estimated in a different population (e.g. the CEPH family parents) from the population from which the disease family is drawn. It is our experience that allele frequencies may vary dramatically between populations. Therefore it may be inappropriate to perform a full likelihood analysis of disease linkage, if the results of that study depend in an essential way on parameters estimated from another population. A somewhat more limited analysis which makes no assumptions concerning allele fre- quencies, such as that given by CRIMAP, may be preferable. If (for reasons of power) a full likelihood analysis is necessary, that analysis should be performed using a range of allele frequencies at each locus, in order to ensure that the results do not depend in an essential way on these frequencies. Finally, it should be noted that whenever some available information is arbitrarily excluded, as is the case with CRI-MAP's analysis of general pedigrees, there is a potential that artifactual biases in parameter estimates can be created. In a preliminary examination of this possibility (cf. Goldgar et al, 1989), we have compared the results of CRI-MAP with those of a full likelihood analysis program (LINKAGE) for all pairs of a set of 20 linked loci, in a data set having 136 pedigrees, most with missing individuals. The information loss (as estimated by comparing the effective number of informative meioses for the two programs) averaged about 15%, and the estimated recombination fractions were almost exactly the same. This suggests that bias is not significant and information loss is small, at least in this particular data set. However, the possibility of bias clearly needs to be considered more thoroughly by means of simulation studies involving multiple loci, and we are in the process of doing this. I would recommend that you check the extent of information loss and/or possible bias in your data set by comparing the results of CRI-MAP to those of a full likelihood program for small sets of loci. GETTING STARTED 1. Put all files on the diskette provided into a single directory on your computer. 2. Compile the source code files (i.e. those whose names end in .c), and link them together with your C stdio and math libraries, producing an executable file called crimap. There are slight differences in the conventions used by some C compilers which may make it necessary for you to modify some of the source files before they will compile correctly; usually these will be evident from the error messages generated by the compiler. For example, if your operating system is case-sensitive (e.g. UNIX), you will need to rename all files in lower case (DOS converts all file names to upper case). We have declared some stdio and math library functions in the crimap functions where they are used; if your compiler objects to this, you will need to delete the offending declarations. Note to VAX C users: there is a bug in the optimizer of the VAX C compiler, ver. 2.3 (and possibly in other versions), which causes it to compile the source file "e_conv.c" incorrectly. When the program is run, this will result in an incorrectly created phase-known data structure in the .dat file (see below). To avoid this, compile e_conv.c with the /NOOPTIMIZE option. All other source files should be compiled with the default /OPTIMIZE option. Also see comments on "underflow" under TECHNICAL NOTES. 3. The command format for running the program is crimap {chromosome number} {option} Example: crimap 7a twopoint The option name must be entered in lower case letters. Chromosome "number" (which may consist of any string of digits, possibly followed by letters -- for example 7p, or 17nf, or 0a) may be replaced by the name of the parameter file (described below); for example, crimap chr7a.par twopoint You must provide a .gen file (named in accordance with the conventions described below, and residing in the same directory as CRI-MAP) which contains the raw genotype data, and run the prepare option first in order to create the other files required by the program. All program output (apart from specific information written to one of the four files described in the next section) is displayed using the "printf" function in C. It will thus be displayed on the terminal (unless redirected to a file by means of commands to the operating system) if the program is run interactively, or written to a log file if the program is run in batch mode. (The latter procedure is the most convenient way to make a copy of the output. In UNIX and some other operating systems, one can simply redirect the output to a file.) The prepare and merge options are the only ones requiring interactive input. 4. As a test run with the chromosome 7 data set, chr7a.gen, provided with the program, use the command crimap 7a prepare to create a .dat file and a .par file for subsequent use by the option all, with the loci 2 8 9 10. (Specify any two of these as the "ordered loci", and the other two as the "inserted loci".) Use default values for the other parameters. NOTE: If the program stops prematurely, displaying the message "Your compiler uses a different size for integers; see documentation for changes that will have to be made in the source code", then you will need to make certain minor modifications in the source files and recompile. See "16 bit integers" under TECHNICAL NOTES, below. If the program stops prematurely, displaying the message "ERROR: ALLOCATION FAILED IN MORECORE", then the program's default memory allocations exceed what your operating system can provide. Reduce the value of nb_our_alloc in the .par file (using prepare, or a text editor) and see "Memory management" under TECHNICAL NOTES, below. Then run CRI-MAP with the command line crimap 7a all The program prints out parameter values, tests all 12 possible orders of these 4 loci, and finally prints out a sorted list of all orders whose log 10 likelihood is at least that of the best order, less 3.0. You should get the answer: 2 10 9 8 -97.937 2 10 8 9 -100.839 The run to evaluate these 12 orders should take less than a minute on most workstations or moderately powerful personal computers (for example, about 50 secs on a Mac II, about 25 secs on a MICROVAX II without a numerical coprocessor, about 8 seconds on a SUN 3/60 or 2 seconds on a 4/60; exact times will depend on the efficiency of the C compiler used as well as the computer). FILE STRUCTURES CRI-MAP uses four different types of files. These all have names of the form chrx.y, where x is the chromosome "number" and the suffix .y is one of .gen the raw genotype data file .par the parameter file .dat the processed data file .ord the orders file, or "orders database" The user must supply the .gen file. The .dat and .par files are created using the option prepare, and must be constructed before using any other CRI-MAP options (except merge). The .ord file is initialized by prepare, but required only for the map building options (build, instant, and quick). You will need to learn about the structures of the .gen and .par files, but can ignore the descriptions of the other file types if you wish. Each file is in ASCII format, and can be edited with a text editor. For readability, the user can insert additional blank or end of line characters into these files if desired, since the C statements which read the files treat any string of such characters as a single "white space" delimiter between data items. All data items are either numbers or character strings containing no blank spaces; they must be separated by white space delimiters. In the following descriptions, brackets {} are used to enclose data items appearing in the file but do not themselves appear in the file. .gen file The data items in this file are as follows: {# of families} {# of loci} {locus name1} {locus name2} . . . (each name consisting of at most 15 characters) For each family: {family ID} (consisting of an arbitrary character string) {# of members} For each member in the family: {ID #} {mother's ID #} {father's ID #} {sex: 0 if female, 1 if male, 3 if unknown} {locus1 allele1}{locus1 allele2}{locus2 allele1}{locus2 allele2} . . . The pedigree structure must be completely specified; this means that for any individual with an ancestor in the pedigree, both parents must be assigned ID #'s and must appear in the file. The family ID may be any string of characters (without embedded blanks) but individual IDs must be numbers. For "original parents" (individuals without ancestors in the pedigree), the mother and father IDs are coded as 0. Alleles must be represented by integers, with missing alleles scored as 0. To handle new mutations in a dominant disease gene correctly, ancestors of the mutant individual should have their disease locus genotypes coded as missing. To handle non-pseudoautosomal X linked loci, a dummy allele number (e.g. 9) must be created and assigned as the second allele for all males in the pedigree. Example: For a data set with two loci, LOCA and LOCB, scored on the single pedigree depicted on the next page, the corresponding .gen file would be 1 2 LOCA LOCB P100 6 1001 0 0 0 1 2 2 2 1002 1001 1009 0 1 2 1 2 1003 0 0 1 1 1 1 2 1004 1002 1003 0 1 2 1 1 1005 1002 1003 0 1 1 2 2 1009 0 0 1 0 0 0 0 Note that the missing maternal grandparent has been assigned an ID of 1009 and included in the file. .par file This file may either be created using the prepare option, or directly using a text editor. The format has been changed effective with version 2.4, and differs somewhat from that of the other files. Each line con- tains a parameter or variable name and its value(s) (separated by spaces), and concludes with an asterisk *. The parameters may appear in any order within the file. If a parameter is omitted (either by omitting the corresponding line altogether, or by omitting the value between the parameter name and the asterisk), then CRI-MAP will automatically assign a default value to the parameter. In the following descriptions, each line ending in an asterisk is displayed as it would appear in the .par file, and gives a parameter name together with its default value(s) (if any). Any number of hap_sys, hap_sys0, and fixed_dist lines may appear in the file. The other parameters should appear only once. (If they appear more than once, only the last specified value will be used). dat_file chrx.dat * gen_file chrx.gen * ord_file chrx.ord * These give the names of the .dat, .gen, and .ord files to be used in the analysis. When the name of the .par file is not of the form chrx.par, the default names for the above files are obtained by replacing the ".par" (which must always appear) in the name of the .par file with ".dat", ".gen", or ".ord", respectively. (This is sometimes useful if one wants to keep all files in a separate directory from the program itself, in which case the name of the .par file passed to crimap may include the full path name; the default names for the other files will then include the same path). nb_our_alloc 3000000 * The initial memory allocation (in bytes). The default value will suffice for most runs, except with a very large number of loci, and in fact is much more than is usually needed for analyzing a small number of loci. If additional memory is needed during the run, it will be allocated automatically (up to the limits of your system); however it is advantageous to choose values for the initial allocation which will suffice for the entire run. See "Memory management" under TECHNICAL NOTES. SEX_EQ 1 * SEX_EQ is 1 if recombination rates are assumed equal in the two sexes, 0 if sex-specific rates are to be allowed in each interval. TOL .01 * The tolerance for determining convergence of the layered EM algorithm; when log 10 likelihoods from successive "phase unknown" iterations increase by less than this amount, iteration terminates. The tolerance used to detect convergence in the "non-informative locus" part of layered EM, and in the option twopoint is TOL/10. In extensive tests using our RFLP CEPH family data sets (several hundred maximum likelihood estimations) the default value .01 was always found to be adequate (it was occasionally not adequate when ordinary EM was used as the search method, instead of layered EM). If you are concerned about the possibility that this may not be stringent enough for for your data set (for example, if the likelihood surface is relatively flat), try using .001 instead; the linear nature of EM convergence guarantees that, if the estimates had not converged with TOL = .01, then a substantial improvement in likelihood should be apparent with the more stringent tolerance. If such an improvement is seen, use successively smaller values for TOL until no further improvement in the likelihood results. PUK_NUM_ORDS_TOL 6 * Applies only to the option build; gives the maximum number of orders allowed in the current map, in the phase unknown part of the analysis. PK_NUM_ORDS_TOL 8 * Similar to PUK_NUM_ORDS_TOL, but applies instead to the phase known analysis in build. If 0, the phase known analysis will be skipped entirely during mapbuilding. PUK_LIKE_TOL 3.0 * The tolerance for discarding locus orders. If the log 10 likelihood of an order is less than the log 10 likelihood of some other order for the same loci by an amount exceeding PUK_LIKE_TOL, that order is discarded (or not printed). Used by mapbuilding options, all, and flipsn. With twopoint, LOD tables are displayed only for locus pairs whose LOD exceeds PUK_LIKE_TOL. PK_LIKE_TOL 3.0 * As above, but applies only to analysis of the phase known data in the option build (and is not used by the other options). use_ord_file 0 * This parameter applies only to the options all and flipsn. When it is 1, the orders generated by those options are prescreened against the .ord file to eliminate orders incompatible with the orders database, prior to computing likelihoods. When use_ord_file is 0, the information in the orders database is not used. write_ord_file 1 * Applies only to the option build. When it is 1, the results of the current build run are used to update the orders database. When it is 0, the orders database will not be updated (but will still be used to prescreen orders during the course of the run). ordered_loci {index # of 1st ordered locus} {index of 2d ordered locus} . . . * inserted_loci {index # of 1st inserted locus} {index of 2d inserted locus} . . . * Note: locus indices start at 0. For all options except twopoint, the "ordered_loci" are assumed to be in their known, unique order; the remaining loci, called the "inserted_loci", are to be placed in the framework defined by the ordered loci. For the options chrompic, fixed, and flipsn, there are no inserted loci. For twopoint, if both ordered_loci and inserted_loci are specified, then LOD tables are only computed for pairs of loci for which one is in the ordered list and the other is in the inserted list. Otherwise the analysis uses all pairs of loci in the specified list (ordered or inserted). hap_sys {index # of 1st locus in system} {index # of 2d locus} . . . * hap_sys0 {index # of 1st locus in system} {index # of 2d locus} . . . * Each haplotyped system is a list of loci (for example, different RFLPs detected by the same probe) which are to be grouped together in an analysis. The first locus in any system is called "primary", and the remaining loci are called "secondary". When the parameter use_haps (see below) is 1, the secondary loci in a system "tag along" with the primary locus whenever ordered sets of loci are constructed. The operations which construct new orders (for example, by inserting a new locus into the map, or by permuting a collection of loci) utilize only the primary loci; once the order is constructed, however, it is "filled out" by inserting secondary loci immediately following the corresponding primary locus, prior to calculating likelihoods and map distances. Secondary loci are automatically deleted from the input lists of ordered loci and inserted loci. Thus any system which is to be included in the analysis must be represented by (at least) its primary locus. For systems specified using hap_sys, the loci within a system are treated as independent in all calculations; i.e. they are not forced to have 0 recombination fraction. In particular, intralocus recombinants between loci in the same system are permitted. For systems specified using hap_sys0, distances within the system are forced to 0 (the program will stop, displaying the message ERROR: 0 likelihood, if there are in fact intralocus recombinants). Example: two haplotyped systems, the first having loci 3, 4, and 5 with distances not forced to 0.0, and the second having loci 9 and 11 with distances forced to 0.0, would be entered in the .par file as hap_sys 3 4 5 * hap_sys0 9 11 * use_haps 1 * When use_haps is set to 1, haplotyping is performed; when it is 0, any haplotyped systems specified in the .par file are ignored (i.e. the input lists of loci are taken as is, and no secondary loci are deleted or inserted). fixed_dist {rec. frac.} {index # of 1st locus} {index # of 2d locus} {sex (optional)} * For the options fixed and chrompic (only), the recombination fraction between a pair of adjacent loci may be held fixed using fixed_dist. If the recombination fraction to be fixed is sex-specific, specify the sex as 0 for female, or 1 for male. If either locus is part of a haplotyped system, it must be the primary locus from its system (it is not possible to force a distance within a haplotyped system, except by using hap_sys0 as described above). Note: any recombination fractions held fixed, either using fixed_dist or hap_sys0, are flagged by an asterisk in the map displayed following the analysis. The last line of the .par file must contain the single word END. Example: if you wished to construct a .par file chr7a.par to perform the all analysis described in GETTING STARTED, above, the following three lines would suffice (since default values are to be used for all other parameters): ordered_loci 2 8 * inserted_loci 9 10 * END .dat file This file, which is created by the program (using the prepare option), not the user, has two data structures, each having the following data items: {# of loci} {# of families} {family id1} {family id2} . . . (may be any character strings not containing blanks) {locus name1} {locus name2} . . . (each locus name consists of at most 15 characters, with no embedded blanks) For each family: {# of chromosomes} {phase chromosomes: each is a character string of length numloci with values 0, 1, or X, denoting the phase} {# of switches} For each switch: {index of locus affected by the switch (in this file, locus indices start at 1, not 0)} {string of length numchroms, with entry 1 if the corresponding chromosome is affected by the switch, 0 otherwise} The first such data structure has the "phase known" data: there is 1 "family" with all the chromosomes in the data set, and 0 switches. All loci of unknown phase are given phase X, as are loci on the chromosomes of children of identical heterozygotes (the latter restriction is necessary to avoid bias in the estimation of recombination fractions, cf. Ott (1985)). The second data structure contains the full phase information for the data set, arranged by families. (For definitions of switches and phase chromosomes, see (Green, 1988)). The "phase known" data structure is used only by the option build. .ord file (For definition of terminology, see description of the build option, below.) {# of orders objects} For each orders object: {# of orders} {# of loci} For each order in the orders object: {index of 1st locus in the order} {index of 2d locus in the order} ... PROGRAM OPTIONS all Finds log 10 likelihoods for all locus orders which result from placing the inserted loci in all possible positions with respect to the ordered loci (see description of the .par file, above, for definition of these terms). The program finds the order with the highest log 10 likelihood L and prints out, in decreasing order of log 10 likelihood, those orders whose associated log 10 likelihoods are greater than or equal to L - PUK_LIKE_TOL. In particular, when the number of ordered loci is two, all possible orders of the specified loci are analyzed. When using this option, keep in mind that the number of orders analyzed is n!/f!, where n is the total number of loci and f is the number of ordered loci (and ! denotes factorial), and that analysis of all 360 orders of a typical set of 6 loci on chromosome 7 takes about 30 mins on a MICROVAX II; thus it is probably impractical to use this option with more than 6 inserted loci, except on a much faster computer. When use_ord_file is set to 1 in the .par file, any orders incompatible with the orders database are eliminated prior to calculating likelihoods, which may make it feasible to use all with larger sets of loci. One convenient use of this option is to position a new locus (which would be specified as the single inserted locus) with respect to a set of loci of known order. build Builds a map by sequential incorporation of loci. This is the main option used for RFLP map construction. To describe the method used in greater detail, define an "orders object" to be a collection of loci, together with a set of permissible orders for those loci. The "orders database" contains a collection of orders objects. At any given stage during the map construction "the map" (designated "current_orders" in the program source code, and output) consists of one such orders object. At the beginning of the run the orders database is empty, and the map consists of the "ordered loci" specified in the .par file. (In the absence of any prior information concerning locus order, such as physical localization data, one usually chooses the "ordered loci" to be a pair of linked and highly informative loci, in order to accelerate the map construction process). The first locus in the list of inserted loci (specified in the .par file) is placed in each possible interval in the map. The resulting locus orders are then tested for compatibility with the database; each order not excluded is subjected to a full maximum likelihood estimation. The order having the highest log 10 likelihood is found, and any order whose log 10 likelihood is less than this one by more than a specified tolerance (PUK_LIKE_TOL or PK_LIKE_TOL) is eliminated. The resulting collection of non-excluded orders form an orders object (called orders_temp in the program output). The database is then updated by (1) appending orders_temp to it; and (2) for each orders object already in the database, eliminating any order which is incompatible with orders_temp. If orders_temp has fewer, or the same number, of orders as does the map, then it becomes the new map and the added locus is deleted from the list of inserted loci. Otherwise the next locus in the list is tried in the same manner. If no locus in the list meets this criterion, the one with the smallest orders_temp is added to the map. Each time a locus is added to the map, the program returns to the beginning of the (revised) list of inserted loci. The orders database is kept in memory during the program run; a copy of it is written to the .ord file every time a new locus is added to the map, so that the program can be stopped and restarted later without loss of the information in the database. Build first analyzes the "phase known" data (see description of .dat file structure), for which maximum likelihood computation is much more rapid than for the full data set. When the map gets too large (i.e. the number of orders in the map exceeds PK_NUM_ORDS_TOL), build proceeds to find a set of uniquely ordered loci, starting with the original ordered loci and adding additional loci which now have a unique placement (during this latter step, the program only makes use of information in the data base; it performs no likelihood calculations). Using this uniquely ordered set of loci as the new "map", the program then proceeds to analyze the full data set, again adding loci sequentially according to the procedure described above. This portion of the program stops when the number of orders in the map exceeds PUK_NUM_ORDS_TOL. Build then again extracts a set of uniquely ordered loci on the basis of information in the orders database, and prints out sex-specific and sex-averaged recombination fractions and the corresponding Kosambi centiMorgans for these loci. For each remaining locus, it prints out the possible placements with respect to the uniquely ordered loci, along with the log 10 likelihoods for each placement. If one starts the map with a pair of unlinked loci, then construction of the map is slower, because in the initial stages the map will consist of two unlinked linkage groups which can be in either orientation with respect to each other so that there are more orders to keep track of (and to place remaining loci in). Conversely, if the initial loci are at 0 recombination fraction (for example, two RFLPs detected by the same probe), then in all subsequent maps the loci will appear in both possible orientations, thereby doubling the number of orders. What is worse, the parts of the program which extract the maximum set of uniquely ordered loci will never get past these two (since no other locus is uniquely placed with respect to them). To avoid this, it is advisable to put these loci later in the list, or (better) to haplotype them using hap_sys or hap_sys0 in the .par file. It is a good idea to exclude relatively uninformative loci from the initial map construction process: they tend to have multiple possible placements, resulting in a large number of orders to test. Their positions can be determined later using all, or instant. With large numbers of loci, it may be necessary to perform several map runs using build until one arrives at "the best" set of uniquely ordered loci. Information in the orders database is cumulative between runs, provided it is not reinitialized using prepare; it is not necessary to reinitialize the orders file when using different sets of loci from the same file. NOTE: The goal of multilocus linkage analysis is to find the locus order having the highest likelihood, and identify alternative orders with comparable likelihoods. Because it is impossible to consider all possible orders for a large set of loci, it is necessary to adopt a strategy which makes decisions on the basis of subsets of the loci. Multiple statistical tests are performed, each with a nonzero probability for Type I error. Because incorrect orders will sometimes have significantly higher likelihoods than the correct one, and the chance of this happening increases with the number of sets of loci which are examined during the map construction process, it is possible that the maximum likelihood order would in fact be rejected by build (or any other published method for map construction, for that matter) because for some subset of the loci, an alternative order has a significantly higher likelihood. In our experience this happens only very rarely (probably because the tests are not independent), and when it does is usually due to errors in the data itself. Nevertheless, three precautions should be taken to guard against the possibility of an erroneous order being adopted by build. First, we usually construct ini- tial maps using a fairly stringent tolerance for rejecting orders -- PUK_LIKE_TOL at least 3.0 -- and then lower it to 2.0 (and omit the phase known part of the analysis) for construction of the final map. Second, it is advisable to check the final order by running flipsn to look at permutations of successive triples or quadruples of loci. Finally, it is a good idea to construct the map several different times, using different starting pairs of loci and a different sequence for adding the loci each time. chrompic Displays the grandparental origins of alleles in each child's chromosomes, for the phase choice having the highest likelihood; gives the relative probability for this phase choice and for the next most likely alternative; gives the number and location of recombinations on each chromosome, and provides, for each chromosome interval defined by the loci, a list of the chromosomes having a recombination in that interval. Flags candidate data errors. The program begins by finding maximum likelihood estimates of the recombination fractions (for the specified locus order); these are then used to find, for each pedigree in the data set, the particular phase choice having the highest likelihood for that pedigree. (The algorithm used for this purpose is described in (Green, 1988)). For each individual with parents in the pedigree, the chromosomes are displayed (maternal first, then paternal). (The individual's ID # appears to the left of the maternal chromosome.) A '0' (resp. '1') means the allele is of known phase and is of grandmaternal (resp. grandpaternal) origin. For alleles of unknown phase, 'o' and 'i' are used instead, denoting the grandparental origin in the most likely phase choice. '-' means that the locus is noninformative. To the right of the chromosome appears the number of recombinations. Beneath the chromosome appear the names of any informative loci which are out of phase with the nearest informative loci on either side; when the loci are closely spaced, these are candidate data errors. Following the chromosome pictures for all the families is a list of recombinant chromosomes. For each pair of loci, the chromosomes are listed in which there is a recombination between those two loci (and no intervening locus is informative). The recombinant chromosomes are designated by family number, individual ID #, and M or P (for maternal or paternal). Following the list of recombinant chromosomes is a list of groups of loci not separated by any crossovers in the data set; and finally, the maximum likelihood recombination fractions which were used in computing probabilities of the phase choices. fixed For a set of loci in a specified order, finds the associated maximum likelihood recombination fractions and map distances and the corresponding log 10 likelihood. The recombination fraction between any pair of adjacent loci may be held fixed by using fixed_dist in the .par file. flipsn For each locus order obtained by permuting an n-tuple (n >= 2) of adjacent loci within an initial (reference) locus order, displays the relative log 10 likelihood (i.e. the log 10 likelihood of the reference order, minus that of the permuted order). For example, flips3 will find relative log 10 likelihoods for all permutations involving triples of adjacent loci. (If n is omitted it is assumed to be 2). When n is 2, relative log 10 likelihoods for all permutations are displayed; for higher values of n, only the permutations whose log 10 likelihoods are at least that of the original order, less PUK_LIKE_TOL, are displayed. Note: the total # of orders which must be evaluated is (m - n + 1)(n! - (n-1)!) + (n - 1)! where m is the total # of loci; thus it is impractical to use this option with large values of n. However, when use_ord_file is set to 1 in the .par file, permutations incompatible with the orders database are eliminated before calculating likelihoods, which will sometimes make it feasible to use larger n's. instant Finds a uniquely ordered set of loci, and the possible placements of the remaining loci with respect to them, using only the information already existing in the .ord file (the .ord file must have been created in a previous build run). Log 10 likelihoods for the different placements of the non-uniquely ordered loci are computed. merge Merges two .gen files, having overlapping sets of families and/or loci. The option is run interactively: in response to prompts, the user enters the names of the two files to be merged, and the name for the merged out- put file. File names must be complete (i.e. the .gen suffix needs to be included). The program compares family IDs and locus names in the two input files, and for families having the same ID compares individual ID #s; it then merges the data from the two files into a single .gen file, taking overlaps into account. A family appearing in both files need not have exactly the same individuals (in which case the output file will have a merged family with all individuals appearing in either input file); the program checks individuals which do occur in both input files to make sure their mother and father ID #s and sex agree, and if they don't displays an error message. If the two files have some loci in common, the genotypes at those loci in the output file for an individual appearing in both input files will be those from the second file. Note: although the command format for using this option is the same as for the other options, the chromosome number is ignored (although it must appear), and no .par file is needed. N.B. Be careful when merging files containing, say, CEPH and disease families that there are no spurious matches between family IDs. For this purpose it is useful to use the fact that family IDs may be arbitrary character strings; so for example CEPH family IDs may be prefixed with "CEPH". prepare This option (which is run interactively) must be run before any of the other options, to create the .dat and .ord files used by them. It will also create a .par file (or modify an existing one), however once you are familiar with the format of the .par file you will probably find it quicker to create and modify that file with a text editor, rather than use prepare repeatedly. Example of usage: crimap 7a prepare Prepare first searches for the .gen and .dat files of the correct name (in the above example, chr7a.gen and chr7a.dat). If there is no .dat file, one is created from the .gen file; if both exist, the program asks whether a new .dat file should be created from the .gen file (select this option if the .gen file is more recent). If neither exists, a message to that effect is given and the program terminates. As the .dat file is being created, family IDs are displayed and any non-inheritances are noted. (These should be checked carefully: a missing data item in the .gen file will generally result in multiple non-inheritances and incorrect family IDs for the subsequent families. If non-inheritances are found, these should be corrected in the .gen file, and prepare run again). Prepare then displays values for several parameter values, including whether to compute sex-averaged or sex-specific likelihoods, and tolerances which control convergence of the maximum likelihood search and the map building process (see description of .par file, above). The displayed values are the defaults (given in the description of the .par file) or the values specified in any preexisting .par file of the same name. The user is then prompted to change any of these values, if desired. Any haplotyped systems, or fixed distances (see description of .par file) in a preexisting .par file are then displayed, and you are prompted to enter any new ones if desired. (prepare cannot be used to modify or delete systems in the preexisting file; you must use a text editor for this). Prepare then displays the names of the other CRI-MAP options, and asks which one it should set up the .par file for. Depending on the option chosen, the user is now prompted to specify the indices of the "ordered loci", and of the "inserted loci". The ordered loci are generally fixed in that order during subsequent analysis, while the inserted loci are tried in all possible positions. When nothing is to be assumed about locus order only two ordered loci should be specified. (Twopoint is a special case; see its description below). For the map-building options (build, instant, and quick), the ordered loci form the initial map. As the default (in this case only), if requested, prepare will sort the loci by informativeness of the phase known data and designate the two most informative loci as the ordered loci, and the remaining loci as inserted. For the options fixed, flipsn, and chrompic, all loci are specified as ordered, and none are inserted. See the description of the .par file, or the relevant option, for further details. If the option build has been specified, prepare initializes an orders file, unless one already exists, in which case it asks whether the existing one should be used. Finally, prepare asks whether to create a new .par file containing the information you have input. You have the option at this point of aborting the creation of a new file, in which case any preexisting .par file will remain intact. quick Similar to instant, except that log 10 likelihoods are not computed. twopoint Performs twopoint linkage analysis for each pair of loci. In the case of sex specific analyses, the output for each pair appears on three lines: the first has the locus names, followed by the maximum likelihood esti- mates, f and m, of the female and male recombination fractions, and the sex-specific LODS ( = log 10 like(f,m) - log 10 like(.5,.5)). The second line has the values of log 10 like(r,m) - log 10 like(.5,.5) for r = .001, .01, .05, and all multiples of .05 up through .5. The last line similarly has the values of log 10 like(f,r) - log 10 like(.5,.5). Only pairs for which the LODS exceed PUK_LIKE_TOL are printed out. (You can get tables for all pairs by setting PUK_LIKE_TOL = 0.0). When SEX_EQ is 1, the sex-averaged recombination fractions, LODS, and log 10 likelihoods are given instead. Haplotyped systems may be used for either locus in a pair. If the loci within a system are not forced to have 0 distance, then maximum likelihood intralocus distances are found simultaneously with the maximum likelihood interlocus distances, and are subsequently held fixed at these values when likelihoods are calculated at the above specified interlocus distances. If both ordered_loci and inserted_loci are specifed in the .par file, tables are only computed for pairs consisting of one locus from each list. This is convenient, for example, when one has just added data for a new set of loci to a file and wishes to look for linkages with the old loci. If only one list (either ordered_loci or inserted_loci) is non-empty, LOD tables are computed for all pairs of loci within that list. TECHNICAL NOTES Memory management Memory is allocated dynamically by the program, as needed; however, since repeated requests of memory from the operating system (via the C stdio library function malloc) are quite timeconsuming, we have adopted a strategy in which two large initial blocks of memory are first allocated. One of these is reserved for the objects comprising the orders database and is apportioned out as needed by the function our_orders_alloc, while the other is used for all other memory requirements and is apportioned by our_alloc. If the amounts initially allocated are insufficient, additional memory is requested from the operating system (via the library function malloc). (When this happens, a statement to that effect is printed). The default initial request of 3Mb (for our_alloc) will suffice for nearly all runs. If the program requests more memory than the operating system can provide, the program terminates with the message "ERROR: ALLOCATION FAILED IN MORECORE". If this happens, try setting the amount of memory specified in the .par file (using the parameter nb_our_alloc) to the maximum amount the operating system will allow you (and see if your system manager can increase that amount by adjusting operating system parameters). If that fails, you may need to reduce the number of loci being analyzed (often, deletion of a single locus having a large number of switches may suffice), or to split a large pedigree into subpedigrees. I would be interested to learn of cases where that is necessary, since it may be possible to improve the switch algebra algorithm to avoid it. Underflow When mapping large numbers of loci, the individual family likelihoods involve products of many factors numerically less than 1, and so can become quite small. If possible, at compilation time you should select an option for representation of floating point numbers which allows them to be as small as possible. With VAX C, this is done by using the /g_float option with the compile command cc, and then linking with the vaxcrtlg library. 16 bit integers C compilers usually represent integer variables with either 16 or 32 bits. The source code provided to you assumes that you have a "32 bit" compiler. If, when you first try to run the program, it terminates prematurely, displaying the message "Your compiler uses a different size for integers; see documentation for changes that will have to be made in the source code", then you have a 16 bit compiler. To get the program to run correctly you must then make the following two changes in the source code and recompile: (1) in the file "defs.h", replace the line typedef int INT; by typedef long INT; (2) Consult the manual for your compiler and determine the name of a library memory allocation function which takes long integers, instead of integers, as arguments (the function malloc takes integer arguments). The name of this function will be something like "mlalloc" (Lightspeed C) or "_halloc" (Microsoft C, ver 4.0 or later). Whatever it is, substitute it for "malloc" throughout the source files our_allo.c and our_orde.c. If your compiler does not have such a function, you will be unable to run CRI-MAP. CHANGES IN VERSION 2.2 - 2.4 The main change in version 2.2 is the incorporation of a more efficient "switch algebra" algorithm which in many cases (particularly for large pedigrees) substantially reduces memory requirements and running time. Thus for many runs, even with a relatively large number of loci, the default memory request of 3 MB (for nb_our_alloc in the .par file) is overly generous. Version 2.3 incorporates a better switch algebra algorithm which further reduces memory and running time. Additional changes have been made which should make this version compatible with most C compilers. Version 2.4 incorporates the following changes: (1) Haplotyping: it is now possible to constrain the loci in a haplotyped system to be at 0 distance from each other for all likelihood calculations. In addition, haplotyping can be used with all of the options; in particular, twopoint LOD tables can be generated for pairs of haplotyped systems. Haplotyped systems are now specified in the .par file (and there is no longer a .hap file). (2) Fixed distances: it is now possible to fix the distance between a pair of adjacent loci, for the purpose of likelihood calculations, with the options fixed and chrompic. (3) The .ord file: it can optionally be used to prescreen locus orders generated by all and flipsn, to reduce the number of orders for which likelihoods need be calculated. (Previously the .ord file was only used by the mapbuilding options). The parameter use_ord_file in the .par file controls the use of the .ord file with these options. (4) twopoint: it is now possible to divide the loci into two different groups, and generate LOD tables only for pairs consisting of one locus from each group (this makes it possible, for example, to conveniently generate LOD tables between "new" loci and "old" loci). The parameter PUK_LIKE_TOL is now used as a LOD cutoff: tables are displayed only for locus pairs whose LOD exceeds PUK_LIKE_TOL. (5) chrompic: chromosomes are now labelled by individual ID # (not position within .gen file); phase designations are now always consistent across pedigrees (not merely across chromosomes) even in the phase unknown case; the relative probabilities for the best phase choice and the next most likely alternative are given; the phase designations now distinguish phase known from phase unknown. (6) merge: a new option, merge, makes it easy to merge two .gen files having overlapping sets of families and/or loci. Other changes in v. 2.4, include: family IDs are now allowed to be arbitrary character strings (not merely integers as in previous versions); the format of the .par file has been changed to make it easier to use; the orders displayed by all are sorted by decreasing log 10 likelihood; output of flipsn is easier to read (the "flipped" loci in each order are highlighted). BIBLIOGRAPHY Barker D, Green P, Knowlton R, Schumm J, Lander E, Oliphant A, Willard H, Akots G, Brown V, Gravius T, Helms C, Nelson C, Parker C, Rediker K, Watt D, Weiffenbach B, Donis-Keller H (1987): Genetic linkage map of human chromosome 7 with 63 DNA markers. Proc. Natl. Acad. Sci. USA 84: 8006-8010. Donis-Keller H, Green P, Helms C, Cartinhour S, Weiffenbach B, Stephens K, Keith T, Bowden D, Smith D, Lander E, Botstein D, Akots G, Rediker K, Gravius T, Brown V, Rising M, Parker C, Powers J, Watt D, Kauffman E, Bricker A, Phipps P, Muller-Kahle H, Fulton T, Ng S, Schumm J, Braman J, Knowlton R, Barker D, Crooks S, Lincoln S, Daly M, Abrahamson J (1987): A genetic linkage map of the human genome. Cell 51, 319-337. Goldgar D, Green P, Parry D, Mulvihill J (1989): Multipoint linkage analysis in Neurofibromatosis Type I: An international collaboration. Am J Human Genetics 44: 6-12. Green P (1988): Rapid construction of multilocus genetic linkage maps. I. Maximum likelihood estimation. (draft manuscript). Kernighan B, Ritchie D (1978): The C programming language. Prentice-Hall. Lander E, Green P (1987): Construction of multilocus genetic linkage maps in humans. PNAS 84, 2363-2367. Ott J (1985): Analysis of human genetic linkage. Johns Hopkins University Press.