**WH** -- A COMPUTER PROGRAM FOR ISOLATION MODEL FITTING--*
Department of Genetics
Nelson Biological Labs
604 Allison Rd.
Piscataway, NJ� 08854-8082
*Some key internal parts of this program derive from the first program for the method that was written by John Wakeley.
This computer program and documentation may be freely copied and used by anyone, provided no fee is charged for it.
WH is a program that fits a simple speciation model, called the Isolation Model, to multilocus DNA sequence data sets. The isolation model assumes the following:
Two species, or populations, from which the data have been sampled, arose from a single ancestral species, or population, some time, t generations, in the past.
The common ancestral species, or population was of constant effective size NA.
The two descendant species (populations) have contant effective sizes N1 and N2, respectively.
Except for the moment of population separation at time t, population sizes are constant.
Since population separation, there is zero gene flow between the species (populations).
All mutations are neutral.
WH implements the methods described in
Wakeley, J., and J. Hey, 1997 Estimating ancestral population parameters. Genetics 145: 847-855.
Wang, R. L., J. Wakeley and J. Hey, 1997 Gene flow and natural selection in the origin of Drosophila pseudoobscura and close relatives. Genetics 147: 1091-1106.
Downloadable Files������������� Return to Contents
�Input File Format��������������� Return to Contents
First Required Value - an integer corresponding to the number of loci.
First Optional Value - The mutation rate per year since divergence of the species, totaled across all of the genes
Second Optional - The number of generations per year for the organisms in the species
- NOTE - the mutation rate per year and the number of generations are not required for the primary analysis. If you desire estimated parameter values in terms of actual population sizes, and years, then estimates of these values should be provided.
- If all that is to be done is a basic fitting of the isolation model, and if no tests are to be done, then each line is required to have 9 items - see DATA LINES below.
- If simulations and statistical tests are to be done, then each line also requires 2 additional items � the population recombination rate estimates for each species
- If linkage disequilibrium tests are also to be done, then there will be an additional 6 items, for a total of 17.
FOR EACH LOCUS, ONE LINE PER LOCUS, IN ORDER:
If Simulations are to be done then then each line should also have:
If tests of Linkage disequilibrium are to be done then each line should also have:
Note if there is not an LD measurement for a locus, -10 is used.
Note on simulations.
The simulation results are sensitive to the amount of recombination. In the published descriptions of these simulations (Wang, Wakeley and Hey, 1997; Kliman et al., 2000) we used the gamma estimator of recombination (Hey and Wakeley, 1997). This estimator tends to have a bias such that the estimates are lower than the expected value of the parameter. The result of having lower recombination is to raise the variance of the observations (of exclusive, shared and fixed variants) and thus to broaden the distribution of test statitics of the fit of the model to the data. In this sense, the tests should be conservative. However this is not guaranteed, and users may want to exam the quality of the fit between the model and their data by considering a range of recombination rates.
Recombination rate estimates are sometime not available for both species. Also they are never available for the common ancestral species. Following is the method of assignment of population recombination rates:
- The program takes 4Nc1i as input (4Nc for species 1, locus i) and then sets
4Nc2i = 4Nc1i theta2/theta1
4NcAi = 4Nc1i thetaA/theta1
-Obtaining 4Nc1i depends on whether one has estimates for species 1 or species 2 or both. If only 4Nc1i is available, then that's it. - If only 4Nc2i is available, then 4Nc1i = 4Nc2i theta1/theta2. If both are available, then 4Nc1i = (4Nc1i + 4Nc2i theta1/theta2)/2.
�Running the Program������������� Return to Contents
The program file should reside either in the same folder as the data file or in a folder automatically searched by the operating system. The program can be run using command line parameters, or by simply typing the name of the program ('wh'). If command line parameters are not used, the program asks for the values of runtime parameters.
The user starts the program simply by going to the folder where the data file exists and typing the name of the program (e.g. 'sites') followed by the� enter key. The program asks several questions about the data file and the� desired analysis. Nearly all commands and options can also be entered using command line parameters.
The program can be started with or without the use of instructions at the command line.
Without command line instructions - simply type �wh� at the prompt.
The program will ask for basic information.
On a PowerPC, clicking on the program icon opens a small window in which command line parameters can be entered. The user can also just hit return at this point and the program will request runtime parameters.
Command Line Parameters:
Type and enter
wh -d'datafilename' -r'resultsfilename' -N'numsims' -L'ldtype' -A'ranseed'
'datafilename' is the name of a data file that is in the same folder as the program
'resultsfilename' is the name of the file that will be created to contain the results. The results file will be generated with an extension of '.wh'
'numsims' is an integer of the the number of simulations to be run. Basic model fitting does not require simulations. If this value is not zero or blank, the data file must have the necessary information on recombination rates.
'ldtype' tells the program what type of linkage disequilibrium value to calculate. This option requires that simulations be done, and that the data file include information on recombination and observed LD values. 'ldtype' must correspond to the measure of linkage disequilibrium used to generate the datafile. 'ldtype' can take on the following values:
'R' invokes the correlation coefficient
'S' invokes r^2, the square of the correlation coefficient
'D' invokes the basic gametic disequilibrium measure D
'P' invokes the D' which is equal to D/ Dmax, where Dmax is the maximum possible value of D given the observed allele frequencies.
'B' invokes the absolute value of D, |D|
'ranseed' is a random number seed and is not generally required unless the user wants to repeat simulations exactly.
Output������������������������� Return to Contents
Output is contained in the results file. There are three main sections: INPUT; MODEL FITTING RESULTS; and SIMULATION RESULTS.
INPUT simply lists in tabular form the data in the input file.
MODEL FITTING RESULTS lists the following:
Basic Parameter Estimates - these are the overall model parameter estimates as described in Wang, Wakeley and Hey (1997)
Expected Polymorphism Levels - this table simply lists the observed numbers of exclusive, shared and fixed polymorphisms together with the expected numbers given the parameter estimates.
Other Calculations - including alternative time parameterizations, and estimates of effective population size if the optional parameters were provided at run time.
SIMULATION RESULTS lists the following:
The results of tests of the quality of fit of the isolation model to the data. Two test statistics are used, the chisquare (as described in Kliman et al., 2000) and the wh statistic which is the original statistics described by Wang Wakeley and Hey (1997).
A table of parameter value estimates, together with simulated 95% confidence intervals and means. The means are useful to assess the bias of the parameter estimates.
A table of observed, and simulated mean levels of exclusive, shared and fixed variants. This table should resemble the table that compares observed and expected polymorphism levels.
Linkage Disequilibrium Analysis - several tables list comparisons between observed levels of disequilibrium among shared polymorphisms, among exclusive polymorphisms, and between shared and exclusive polymorphisms.
�Program Limitations������������ Return to Contents
For simulations, the program can only handle a total sample size for each locus of 32. If the program is compiled under Microsoft Visual C++ (as the distributed Win32 version is) then it can makes use of a compiler extension and can handle total per locus sample sizes of 64.
For basic model fitting, without simulations, larger sample sizes can be used.
During simulations, recombination within a locus can occur only between sequence segments. The program has been compiled with 50 segments per sequence, which should be more than sufficient for most data sets. However it is possible that this will not be sufficient for loci with long sequences and high amounts of recombination.
�Literature Cited��������������� Return to Contents
Hey, J., and J. Wakeley, 1997 A coalescent estimator of the population recombination rate. Genetics 145: 833-846.
Kliman, R. M., P. Andolfatto, J. A. Coyne, F. Depaulis, M. Kreitman et al., 2000 The population genetics of the origin and divergence of the Drosophila simulans complex species. Genetics 156: 1913-31.
Wakeley, J. and J. Hey. 1997 Estimating ancestral population parameters.Genetics 145, 847-855.
Wang, R. L., J. Wakeley and J. Hey, 1997 Gene flow and natural selection in the origin of Drosophila pseudoobscura and close relatives. Genetics 147: 1091-106.
This page was last changed September 04, 2013