Copyright  2003 by the Genetics Society of America

Modeling Linkage Disequilibrium and Identifying Recombination Hotspots
Using Single-Nucleotide Polymorphism Data
Na Li* and Matthew Stephens†,1*Department of Biostatistics and †Department of Statistics, University of Washington, Seattle, Washington 98195
Manuscript received January 30, 2003
Accepted for publication August 11, 2003
ABSTRACT
We introduce a new statistical model for patterns of linkage disequilibrium (LD) among multiple SNPs
in a population sample. The model overcomes limitations of existing approaches to understanding,summarizing, and interpreting LD by (i) relating patterns of LD directly to the underlying recombination
process; (ii) considering all loci simultaneously, rather than pairwise; (iii) avoiding the assumption thatLD necessarily has a “block-like” structure; and (iv) being computationally tractable for huge genomic
regions (up to complete chromosomes). We examine in detail one natural application of the model:
estimation of underlying recombination rates from population data. Using simulation, we show that in
the case where recombination is assumed constant across the region of interest, recombination rateestimates based on our model are competitive with the very best of current available methods. More
importantly, we demonstrate, on real and simulated data, the potential of the model to help identify andquantify fine-scale variation in recombination rate from population data. We also outline how the model
could be useful in other contexts, such as in the development of more efficient haplotype-based methods
for LD mapping.

L

INKAGE disequilibrium (LD) is the nonindependence, at a population level, of the alleles carriedat different positions in the genome. The patterns of
LD observed in natural populations are the result of
a complex interplay between genetic factors and the
population’s demographic history. In particular, recombination plays a key role in shaping patterns of LD in a
population. When a recombination occurs between twoloci, it tends to reduce the dependence between the alleles
carried at those loci and thus reduce LD. Although
recombination events in a single meiosis are relatively
rare over small regions, the large total number of meiosesthat occurs each generation in a population has a substantial cumulative effect on patterns of LD, and somolecular data from population samples contain valuable information on fine-scale variations in recombination rate.
Despite the undoubted importance of understanding
patterns of LD across the genome, most obviously because of its potential impact on the design and analysisof studies to map disease genes in humans, most current
methods for interpreting and analyzing patterns of LD
suffer from at least one of the following limitations:
1. They are based on computing some measure of LD
defined only for pairs of sites, rather than considering
all sites simultaneously.

1
Corresponding author: Department of Statistics, University of Washington, Padelford B-313, Stevens Way, Seattle, WA 98195.
E-mail: stephens@stat.washington.edu

Genetics 165: 2213–2233 (December 2003)

2. They assume a “block-like” structure for patterns ofLD, which may not be appropriate at all loci.
3. They do not directly relate patterns of LD to biological mechanisms of interest, such as the underlying
recombination rate.
As an example of the limitations of current methods,
consider Figure 1, which graphically shows pairwise LD
measures for six simulated data sets, simulated undervarious models for heterogeneity in the underlying recombination rate. The reader is invited to speculate
on what the underlying models are in each case—theanswer appears in the Figure 8 legend. In each of the six
figures one can identify by eye, or by some quantitative
criteria (e.g., Daly et al. 2001; Olivier et al. 2001; Wanget al. 2002), “blocks” of sites, such that LD tends to behigh among markers within a block. In some cases there
might also be little LD between markers in different
blocks, which might be interpreted as evidence for variation in local recombination rates: low recombination
rates within the blocks and higher rates between theblocks. Indeed, Jeffreys et al. (2001) have shown, using
sperm typing, that in the class II region of the MHC,
variations in local recombination rate are indeed responsible for block-like patterns of LD. However, without this type of experimental confirmation, which iscurrently technically challenging and time consuming,
it is difficult to distinguish between blocks that arise
due to recombination rate heterogeneity and blocks
that arise due to chance, perhaps through chance clustering of recombination events in the ancestry of theparticular sample being considered (Wang et al. 2002).

2214

N. Li and M. StephensFigure 1.—Plots of LD measurement, |D⬘|, (bottom right diagonal)
and P value for Fisher’s exact test
(top left diagonal) for every pairof sites with minor allele frequency ⬎0.15 in data sets simulated under varying assumptions
about variation in the local recombination rate. Details of the models used to simulate each data set
appear in the Figure 8 legend,
which is based on the same six data
sets.

The ability to distinguish between these cases would ofcourse be interesting from a basic science standpoint—
for example, in helping to identify sequence characteristics associated with recombination hotspots. In addition,
it would have important implications for the design andanalysis of LD mapping studies. For example, it would
help in predicting patterns of variation at sites that havenot been genotyped (perhaps sites influencing susceptibility to a disease), and it would provide some indicationof whether block structures observed in one sample
are likely to be replicated in other samples—a crucial
requirement for being able to select representative “tag”
single-nucleotide polymorphisms (SNPs; Johnson et al.
2001) on the basis of LD patterns observed in somereference sample.
In this article we introduce a statistical model for LD
that overcomes the limitations of existing approaches
by relating genetic variation in a population sample tothe underlying recombination rate. We examine in detail one natural application of the model: estimation ofunderlying recombination rates from population data.
Using simulation, we show that in the case in which
recombination is assumed constant across the regionof interest, recombination rate estimates based on ourmodel are competitive with the very best of current
available methods. More importantly, we demonstrate,
on real and simulated data, the potential of the modelto help identify and quantify fine-scale variation in recombination rate (including “recombination hotspots”)
from population data.
Although we focus here on estimating recombination
rates, we view the model as being useful more broadly
in interpreting and analyzing patterns of LD across multiple loci. In particular, as we outline in our discussion,
the model could be helpful in the development of more
efficient haplotype-based methods for LD mapping, along

the lines of, for example, McPeek and Strahs (1999),
Morris et al. (2000), and Liu et al. (2001).
MODELS

Background: The most successful current approaches
to constructing statistical models relating genetic variation to the underlying recombination rate (and to othergenetic and demographic factors) are based on the coalescent (Kingman 1982) and its generalization to include recombination (Hudson 1983). Although theseapproaches are based on rather simplistic assumptions
about the demographic history of the population fromwhich individuals were sampled and about the evolutionary processes acting on the genetic region being
studied, they have nonetheless proven useful in a variety
of applications. In particular, they provide a helpful simulation tool (e.g., software described in Hudson 2002),
allowing more realistic data to be generated under various assumptions about underlying biology and demography, and hence aid exploration of what patterns of
LD might be expected under different scenarios (Kruglyak 1999; Pritchard and Przeworski 2001).
Despite the ease with which coalescent models can be
simulated from, using these models for inference remains
extremely challenging. For example, consider the problem of estimating the underlying recombination rate ina region, using data from a random population sample.
It follows from coalescent theory that population samples contain information on the value of the productof the recombination rate c and the effective (diploid)
population size N, but not on c and N separately. It has
therefore become standard to attempt to estimate thecompound parameter ␳ ⫽ 4Nc, and several methods
have been proposed. Some (e.g., Griffiths and Marjoram 1996; Kuhner et al. 2000; Nielsen 2000; Fearn-

Linkage Disequilibrium and Recombination

head and Donnelly 2001) try to make use of all themolecular data available. However, although such methods have been applied successfully to small regions andnonrecombining parts of the genome (Harding et al.
1997; Hammer et al. 1998; Kuhner et al. 2000; Nielsen
2000; Fearnhead and Donnelly 2001), for even moderate-sized autosomal regions (e.g., a few kilobases in humans) they become computationally impractical (Fearnhead and Donnelly 2001). Other methods, many ofwhich are considered by Wall (2000), make use of only
summaries of the data, substantially reducing computational requirements at the expense of some loss in efficiency.
More recently, Hudson (2001) and Fearnhead andDonnelly (2002) proposed “composite-likelihood” methods for estimating ␳ over moderate to large genomic
regions. Hudson’s method is based on multiplying together likelihoods for every pair of sites genotyped, inwhich these pairwise likelihoods are computed via simulation, assuming an “infinite-sites” mutation model (i.e.,
no repeat mutation). This method has been modified
by McVean et al. (2002) to allow for repeat mutation.
Fearnhead and Donnelly’s method is based on dividing
data on a large region into smaller regions and multiplying likelihoods obtained for each smaller region.
These methods, together with the best of the summarystatistic-based methods of Wall (2000), appear to be
the most accurate of existing methods for estimating
recombination rates from patterns of LD over moderate
to large genomic regions. None of these methods, as
currently implemented, allows explicitly for variation inrecombination rate along the region under study.
A new model: Here we describe a new model for LD,
which enjoys many of the advantages of coalescent-based
methods (e.g., it directly relates LD patterns to the underlying recombination rate) while remaining computationally tractable for huge genomic regions, up to entire
chromosomes. Our model relates the distribution ofsampled haplotypes to the underlying recombination
rate, by exploiting the identity
Pr(h 1, . . . , h n |␳) ⫽ Pr(h 1|␳) Pr(h 2|h 1; ␳) . . . Pr(h n|h 1, . . . , h n⫺1; ␳),(1)where h 1, . . . , h n denote the n sampled haplotypes, and␳ denotes the recombination parameter (which may
be a vector of parameters if the recombination rate isallowed to vary along the region). This identity expresses
the unknown probability distribution on the left as a product of conditional distributions on the right. For simplicity
we often use the notation ␲ to denote these conditional
distributions. While the conditional distributions are notcomputationally tractable for models of interest, they
are amenable to approximation, as we describe below.
Our strategy is to substitute an approximation for these
conditional distributions (␲ˆ , say) into the right-handside of (1), to obtain an approximation to the distribution of the haplotypes h given ␳:

2215

Pr(h 1, . . . , h n|␳) ⬇ ␲ˆ (h 1|␳)␲
ˆ (h 2|h 1; ␳) . . . ␲ˆ (h n|h 1, . . . , h n⫺1; ␳).

(2)
We refer to this as a “product of approximate conditionals” (PAC) model and to the corresponding likelihood
as a PAC likelihood, which we denote LPAC. Explicitly,
LPAC(␳) ⫽ ␲ˆ (h 1|␳)␲
ˆ (h 2|h 1; ␳) . . . ␲ˆ (h n|h 1, . . . , h n ⫺1; ␳).

(3)
Similarly, we refer to the value of ␳ that maximizes LPAC
as a maximum PAC likelihood estimate for ␳ and denote itby ␳ˆ PAC.
The utility of the model (3) will naturally depend on theuse of an appropriate approximation for the conditional
distribution ␲. This approximation should be designed
to answer the following question: if, at a particular locus,
in a random sample of k chromosomes from a population, we observe genetic types h 1, . . . , h k , what is theconditional distribution of the type of the next sampled
chromosome, Pr(hk⫹1|h1, . . . , hk)? We are aware of threeforms for ␲ in the literature, each of which attempts toanswer this question under different assumptions forthe genetic model underlying the loci being studied.
The first and best known comes from the Ewens sampling
formula (Ewens 1972). This arises from considering aneutral locus in a randomly mating population, evolving
with constant (diploid) size N and mutation rate ␮ per
generation, and assuming an “infinite-alleles” mutation
model, in which each mutation creates a novel (previously unseen) haplotype. Under these idealized conditions, if we let ␪ ⫽ 4N␮, then with probability k/(k ⫹␪) the k ⫹ 1st haplotype is an exact copy of one of thefirst k haplotypes chosen at random; otherwise it is a
novel haplotype. Although the assumptions underlying
this formula will never hold in practice, it does capture
the following properties that we would expect to holdmore generally:i. The next haplotype is more likely to match a haplotype that has already been observed many times
rather than one that has been observed less frequently.
ii. The probability of seeing a novel haplotype decreases as k increases.
iii. The probability of seeing a novel haplotype increases as ␪ increases.
However, for modern molecular data, and for sequence data and SNP data in particular, it fails to capture the following two properties:
iv. If the next haplotype is not exactly the same as anexisting (i.e., previously seen) haplotype, it will tend
to differ by a small number of mutations from anexisting haplotype, rather than to be completely
different from all existing haplotypes.
v. Due to recombination, the next haplotype will tend
to look somewhat similar to existing haplotypes overcontiguous genomic regions, the average physical2216

N. Li and M. StephensFigure 2.—Illustration of how ␲A(hk⫹1|h1, . . . , hk) builds
hk⫹1 as an imperfect mosaic of h1, . . . , hk. This illustrates thecase k ⫽ 3 and shows two possible values (h4A and h4B) for h4,
given h1,h2,h3. Each of the possible h4’s can be thought of as
having been created by “copying” (imperfectly) parts of h1,h2,
and h3. The shading in each case shows which haplotype wascopied at each position along the chromosome. Intuitively wethink of h4 as having recent shared ancestry with the haplotype
that it copied in each segment. We assume that the copying
process is Markov along the chromosome, with jumps (i.e.,changes in the shading) occurring at rate ␳/k per physicaldistance. Thus the more frequent jumps in h4B suggest a higher
value of ␳ than do the less frequent jumps in h4A. Note thatfor very large values of ␳ the loci become independent, as
they should. Each column of circles represents a SNP locus,
with black and white representing the two alleles. The imperfect nature of the copying process is exemplified at the third
locus, where h4A and h4B have the black allele, although they
copied h2, which has the white allele. In practice, of course,
the shading is not observed, and so to compute the probabilityof observing a particular h4 we must sum over all possible
shadings. The Markov assumption allows us to do this efficiently, using standard methods for hidden Markov models,
as described in appendix a.

length of these regions being larger in areas of the
genome where the local rate of recombination islow.
Stephens and Donnelly (2000) suggested a form for␲ that captures properties i–iv above. In their suggested
form for ␲, the next haplotype differs by M mutations
from a randomly chosen existing haplotype, where M hasa geometric distribution with Pr(M ⫽ 0) ⫽ k/(k ⫹ ␪) (sothat it reproduces the Ewens sampling formula in thespecial case of the infinite-alleles mutation model).
Thus the next haplotype is a (possibly imperfect) “copy”
of a randomly chosen existing haplotype.
Fearnhead and Donnelly (2001; henceforth FD)
extended this form for ␲ to also capture property v
above. In FD’s approximation, the k ⫹ 1st haplotype ismade up of an imperfect mosaic of the first k haplotypes,
with the size of the mosaic fragments being smaller for
higher values of the recombination rate.
Here we use two new forms for ␲ that also capture
properties i–v above. The first, described in detail inappendix a and illustrated in Figure 2, which we denote␲A, is a simplification of FD’s approximation that is easierto understand and slightly quicker to compute. [N. Patterson (personal communication) has independently

suggested a similar simplification.] The second, which we
describe in detail in appendix b and denote ␲B, is a slight
modification of ␲A, developed using empirical results from
Figure 3 to produce a likelihood LPAC that gives moreaccurate estimates of ␳. Where necessary, we denote the
PAC likelihoods and maximum PAC-likelihood estimatescorresponding to ␲A (respectively, ␲B) by LPAC-A and ␳ˆPAC-A(respectively, LPAC-B and ␳ˆPAC-B).
A key property of both ␲A and ␲B is that they areeasy and fast to compute. Unlike the Ewens sampling
formula, but like the approximations of Stephens andDonnelly (2000) and FD, neither corresponds exactly
to the actual conditional distribution under explicit assumptions about population demography and the evolutionary forces on the locus under consideration. Indeed,no closed-form expressions for ␲, based on such explicit
assumptions, and capturing iv or v, are known. However,the suggested forms for ␲ were motivated by consideringboth the Ewens sampling formula and the underlying
genealogy (or, in the case with recombination, genealogies) relating a random sample of haplotypes from aneutrally evolving, constant-sized panmictic population.
As such, it may be helpful to view them as approximations to the (unknown) true conditional distribution
under these assumptions. In particular, there are certain
aspects of many real populations (e.g., population
expansion or population structure) and biological factors (e.g., gene conversion and selection) that these
forms for ␲ do not attempt to capture. For some applications this may not matter very much. For others it may
be necessary to develop forms for ␲ that do capture
these aspects—a point we return to in the discussion.
An unwelcome feature of the PAC likelihoods corresponding to our choices of ␲—and indeed the forms
for ␲ from Stephens and Donnelly (2000) and FD—isthat they depend on the order in which the haplotypes
are considered. In other words, although these likelihoods each correspond to a valid probability distribution on the haplotypes, these probability distributions
do not enjoy the property of exchangeability that wewould expect to be satisfied by the true (unknown)
distribution. Practical experience, and theory in Stephens and Donnelly (2000; their Proposition 1, partd), suggests that this problem cannot be rectified bymaking a simple modification to ␲. Although in principle the dependence on ordering could be removed byaveraging the PAC likelihood over all possible orderingsof the haplotypes, in practice this would require a sum
over n! terms, which is infeasible even for rather smallvalues of n. Instead, as a pragmatic alternative solution,
we propose to average LPAC over several random ordersof the haplotypes. Unless otherwise stated, all results
reported here were obtained by averaging over 20 random orders. In our experience, the performance ofthe method is not especially sensitive to the numberof random orders used—results based on 100 random
orders gave qualitatively similar results, and results

Linkage Disequilibrium and Recombination

2217

Figure 3.—Histograms of the error Err(␳, ␳ˆ PAC-A) ⫽ log10(␳ˆ PAC-A/␳), based on 100 data sets simulated from the standard coalescentmodel with n ⫽ 50 haplotypes and S ⫽ 50 segregating sites. The values of ␳ are (a) ␳ ⫽ 5, (b) ␳ ⫽ 25, and (c) ␳ ⫽ 500.
Superimposed curves are normal densities with the same mean and standard deviation as the 100 values making up the histogram.
These results, as well as those in Figure 4 and Table 1, are based on averaging the likelihoods over 10 random orders of the
haplotypes.

based on a single random order were often not muchworse (data not shown). It is, however, important thatwhen comparing likelihoods for different values of ␳,
the same set of random orders should be used for eachvalue of ␳.

ESTIMATING CONSTANT RECOMBINATION RATEIn this section we consider estimating the recombination rate when it is assumed to be constant across theregion of interest. More precisely, we assume that crossovers in a single meiosis occur as a Poisson process ofconstant rate c per unit (physical) distance and consider estimating the scalar parameter ␳ ⫽ 4Nc. We firstuse simulated data to examine the properties of theestimator ␳ˆ PAC-A, corresponding to the conditional distribution ␲A described in appendix a, under what we callthe “standard coalescent model”: a constant-sized, panmictic population with an infinite-sites mutation model.
We show that, although quite accurate, ␳ˆPAC-A exhibits asystematic bias. We use the empirical results to develop a
modified conditional distribution, ␲B (described in detailin appendix b), whose corresponding estimator, ␳ˆPAC-B,
exhibits considerably less bias and is more accurate. We
compare the performance of models based on both ␲Aand ␲B with results from other methods.
Properties of the point estimate ␳ˆ PAC: We used theprogram mksample (Hudson 2002) to simulate data
sets consisting of samples of SNP haplotypes from thestandard coalescent model for various values of:
1. The number n of haplotypes in the sample
2. The number S of markers typed
3. The value of ␳ (we measure physical distance so thatthe total physical length of each simulated haplotype

equals 1.0. Thus the value of ␳ is also the total valueof ␳ across the region).
For each data set we found ␳ˆ PAC-A by numerically maximizing the PAC likelihood (using a golden bisection
search method; Press et al. 1992) and compared it withthe true value of ␳ used to generate the data.
It seems natural to measure the error in estimatesfor ␳ on a relative, rather than an absolute, scale. For
example, Wall (2000) reported the frequency withwhich different methods for estimating ␳ gave estimates
within a factor of 2 of the true value, and both FD and
Hudson (2001) examine the distribution of the ratio␳ˆ /␳ for their estimates of ␳ˆ and the deviation of thisratio from the “optimal” value of 1. A problem withworking with this ratio directly is that it tends to penalize
overestimation more heavily than underestimation. For
example, overestimating ␳ by a factor of 10 gives a larger
deviation from 1 than underestimating ␳ by a factor of
10. To avoid this problem, we quantify the relative errorof an estimate ␳ˆ for ␳ by Err(␳, ␳ˆ ) ⫽ log10(␳ˆ /␳). This
gives, for example, an error of 0 if ␳ˆ ⫽ ␳, an error of 1
if ␳ˆ overestimates ␳ by a factor of 10, and an error of⫺1 if ␳ˆ underestimates ␳ by a factor of 10.
We note that Err(␳, ␳ˆ ) can also be viewed as the error(on an absolute scale) in estimating log10(␳) by log 10(␳ˆ ).
Thus, if the usual asymptotic theory for maximum-likelihood estimation applies for estimation of log10(␳) in
this setting (which, as discussed in FD, it may not), thenfor the actual maximum-likelihood estimate (MLE) ␳ˆ MLE
of ␳, Err(␳, ␳ˆ MLE) would be normally distributed asymptotically, centered on 0. Optimistically, we might therefore hope that for sufficiently large data sets (large interms of the number of haplotypes, the number of markers, or both) Err(␳, ␳ˆ PAC-A) might be approximately normally distributed, centered on 0. In our simulations, we2218

N. Li and M. Stephensfound that for some combinations of n, S, and ␳ thisdid indeed appear to be the case (e.g., Figure 3b), butthat for other combinations, although the distribution
often appeared close to normal, it was centered aroundsome nonzero value (e.g., Figure 3, a and c), indicatinga systematic tendency for ␳ˆ PAC-A to over- or underestimate
␳. We refer to the median of Err as the “bias” [of log(␳ˆ PAC-A) in estimating log(␳)]. Although bias is usually
defined as a mean error, this is not particularly helpful
here since the mean is often heavily influenced by asmall number of very large values and may even be
infinite in some cases (see also FD). We therefore follow
previous authors, including Hudson (2001) and FD, inconcentrating on the behavior of the median, rather
than the mean, of the error.
Despite the biases evident in Figure 3, a and c, ␳ˆ PAC-A
gives reasonably accurate estimates of ␳. For example,
even in Figure 3c, which shows one of the most extreme
biases that we observed in our simulations, the bias
corresponds to underestimating ␳ by approximately afactor of 2, and ␳ˆ PAC-A is within a factor of 2 of the true
value of ␳ in 68% of cases. Although in many statisticalapplications estimates within a factor of 2 of the truth
would not be considered particularly helpful or impressive, in this setting this kind of accuracy is often noteasy to achieve (see, for example, Wall 2000).
We performed extensive simulations to better characterize the bias noted above and found that, although
the bias depends on all three variables (n, S, and ␳), it
is especially dependent on the average spacing betweensites. More specifically, for fixed n and S we observed astriking linear relationship between the bias and the logof the average marker spacing (Figure 4). This linear
relationship was also apparent for data simulated under
an assumption of population expansion (data notshown). The slope of the linear relationship is negativein each case, indicating a tendency for ␳ˆ PAC-A to overestimate ␳ when the markers are very closely spaced andunderestimate ␳ when the markers are far apart. As thenumber of sampled haplotypes increases, both the slope
and intercept of the line appear to get closer to 0 (Table
1). On the basis of these empirical results we can modify
␲A to reduce the bias of the point estimates (see appendix b for details). The improved performance of thismodified conditional distribution, which we denote ␲B,
is illustrated in the next section.
Figure 4 also illustrates the effect of varying parametervalues on the variability of point estimates. As might
be expected, the variance of the error reduces withincreased sample size and increased number of sites,with the latter providing the more substantial decrease.
For example, doubling the number of sites from 50 to
100 roughly halved the variance of the error in mostcases, while doubling the number of individuals from50 to 100 resulted in much smaller decreases. For a
fixed sample size and number of sites, the variance ofthe error decreases as the spacing between sites grows.

This may be due to the fact that for larger spacings
more recombination events occur, increasing the relative accuracy with which ␳ can be estimated, although
we would not expect this pattern to continue indefinitely
as the marker spacing is increased beyond the range
considered here.
Comparison of point estimates with other methods:
Hudson (2001) introduced a composite-likelihood method
for estimating ␳, based on multiplying together the likelihood computed for every pair of SNPs. He compared
the performance of this method with others in the literature (Hudson and Kaplan 1985; Hudson 1987; Hey
and Wakeley 1997; Wall 2000), under the standard
coalescent model, and found it to be as good as, or
better than, the best of these. We compared the results
reported by Hudson (2001) for his maximum composite-likelihood estimate, ␳ˆ CL, with the results for ␳ˆ PAC-A and␳ˆ PAC-B on data sets simulated under the same conditions(Figure 5). For data sets with small numbers of SNPs(ⱕ ⵑ12) ␳ˆ CL provides the most accurate estimates of ␳,
although all three methods struggle to produce reliable
estimates. For larger numbers of SNPs both ␳ˆ PAC-A and␳ˆ PAC-B tend, desirably, to exhibit less variability than ␳ˆ CL.
Further, ␳ˆ PAC-B exhibits little or none of the bias present
in ␳ˆ PAC-A and provides the most accurate estimates of ␳.
The superior performance of the pairwise compositelikelihood method for data sets with small numbers of
SNPs is perhaps not surprising—indeed, for data setswith only two SNPs ␳ˆ CL is precisely the maximum-likelihood estimate for ␳. However, we note that almost allof the improvement in accuracy comes from the increase in the 10th percentile of the estimator toward
the true value, rather than from a decrease in the 90th
percentile. One possible explanation for this is that ␳ˆ CL
uses a likelihood based on an infinite-sites mutation
model (i.e., assumes no repeat mutation) and so is able
essentially to rule out very small values for ␳ if there iseven one pair of sites at which all four gametes arepresent. (The effect of this may be compounded by thefact that ␳ˆ CL was found by maximizing over a grid of
possible values, which forces all nonzero estimates of ␳
to be above some threshold.) Our estimator does notmake the infinite-sites assumption and so will be more
inclined to estimate very small values of ␳, possibly leading to occasional substantial underestimates. Since inreal data it will typically be unclear whether or not theinfinite-sites assumption holds, the advantage of ␳ˆ CL foreven small numbers of sites is perhaps less clear-cut than
it appears in Figure 5.
We used the same simulated data to examine theaccuracy of estimates of ␳ obtained by the methods
described by Kuhner et al. (2000) and FD, both of which
use computationally intensive Monte Carlo proceduresto attempt to approximate the full coalescent likelihood.
The computational complexity of these approaches increases with what might be called “the total value of ␳
across the region,” or “per-locus ␳,” which we denote ␳˜Linkage Disequilibrium and Recombination

2219

Figure 4.—Box plots showing the relationship of the bias to the average marker spacing. For each combination of parameters,
100 data sets each were simulated under the standard coalescent model. The parameters involved are: the number of haplotypes
in each sample, n ⫽ 20,50,100,200; the number of segregating sites, S ⫽ 20,50,100; and the average marker spacing, ␳/S ⫽0.1,0.5,1.0,5.0, and 10.0. In humans a marker spacing of ␳/S ⫽ 0.5 corresponds to ⵑ1 kb between markers. The unlabeled tickmarks on the y-axis correspond to ␳ˆ PAC-A ⫽ ⫾2␳.

(more precisely, in our notation ␳˜ ⫽ ␳L, where L is thephysical length of the region). Results from FD suggestthat even for small values of ␳˜ (⬍3, say), the approximate-likelihood curves obtained by these methods may
be poor approximations to the actual likelihood curve,
and so it seems unlikely that the curves will be accurate

TABLE 1
The intercepts and slopes of the linear relationshipbetween log10(␳ˆ PAC/␳) and log10(spacing)
Intercept

Slope

n/S

20

50

100

20

50

100

20
50
100
200

⫺0.16
⫺0.12
⫺0.12
⫺0.10

⫺0.12
⫺0.07
⫺0.06
⫺0.05

⫺0.09
⫺0.06
⫺0.04
⫺0.02

⫺0.18
⫺0.16
⫺0.09
⫺0.06

⫺0.21
⫺0.21
⫺0.17
⫺0.14

⫺0.26
⫺0.24
⫺0.21
⫺0.17

See also Figure 4.

for larger ␳˜ . However, point estimates based on thesemethods could still be accurate, if the maximum of theapproximate-likelihood curve occurs in about the rightplace. To investigate this possibility, we applied bothmethods, using ⵑ1 day of CPU time per method per
data set (compared with ⵑ30 sec per data set for ␳ˆPAC-B), to
10 of the data sets simulated with ␳˜ ⫽ 40. Computational
considerations make a more comprehensive simulation
study inconvenient. Each of the methods was run with
␪ fixed at the value used to simulate the data, giving
them some advantage over how they could be used inpractice. Nevertheless, neither method produced point
estimates of ␳˜ as accurate as those from ␳ˆ PAC-B (Table 2).
Of the two full-likelihood schemes, the maximum of
the likelihood curve obtained by infs was consistently
closer to the true value of ␳˜ than was the maximum of
the likelihood curve obtained by Recombine. Indeed,
the estimates obtained from Recombine were often
close to an order of magnitude smaller than the true
value of ␳˜ , which raises a danger that when the method
is applied to real data (for which the value of ␳˜ is of

2220

N. Li and M. StephensFigure 5.—Comparisonof ␳ˆ PAC-A and ␳ˆ PAC-B with Hudson’s pairwise compositelikelihood estimator ␳ˆ CL(Hudson 2001) on data setsof n ⫽ 50 haplotypes simulated from the neutral infinite-sites model. The data sets
were simulated with haplotypes of physical length L ⫽4, 8, 12, 16, 20, 24, 28, 32,
36, 40, 44, and 48 units, with␳ ⫽ 1/unit physical length,and ␪ ⫽ 1⁄4/unit physicallength. (With these parameters the expected numberof SNPs in each data set isapproximately equal to thephysical length of the haplotypes.) The results for ␳ˆCL come from Hudson (2001) and were kindly provided by R. Hudson. The resultsfor ␳ˆPAC-A and ␳ˆPAC-B are based on 1000 data sets that we simulated for each set of parameters using the program mksample (Hudson
2002). (We discarded the few simulated data sets that had only one SNP.) The panels are (a) ␳ˆCL, (b) ␳ˆPAC-A, and (c) ␳ˆ PAC-B. In eachpart, the solid line is the median of Err(␳ˆ , ␳) ⫽ log10(␳ˆ /␳) and the dashed lines are the 10 and 90% quantiles.

course not known) the user might be misled into thinking that the value of ␳˜ is small enough for the method
to produce reliable results. Results from longer runs ofinfs, taking ⵑ5 days of CPU time each, produced improved
results, competitive with ␳ˆPAC-B (data not shown).
Properties of PAC likelihood curves: Construction ofconfidence intervals: We examined the coverage properties of confidence intervals (C.I.’s) constructed fromthe PAC-likelihood curve in two ways:

TABLE 2Comparison of ␳ˆ PAC-B with estimates of ␳˜ from Recombine andinfs, for 10 data sets simulated with ␳˜ ⫽ 40, ␪ ⫽ 10, n ⫽ 50
Data set
1
2
3
4
5
6
7
8
9
10
MedianMedian |Err(␳ˆ , ␳)|

Recombineinfs

␳ˆ PAC-B

13
9
24
10
10
5
4
8
7
12
10
0.62

26
21
27
29
27
33
25
23
53
27
27
0.17

57
24
30
26
34
42
45
38
50
29
36
0.11

Both infs and Recombine were run with ␪ fixed at its true
value. infs was run for 20,000 iterations with five driving valuesfor ␳˜ (10, 30, 40, 50, and 60). The effective sample size (ESS)at the MLE is always ⬍4, indicating that infs had very little
confidence in its estimated likelihood curve (and the estimated 95% C.I.’s failed to include the true ␳˜ in all but onecase). Recombine was run with five short runs of 20,000 iterations and one long run of 1 million iterations, using threeheating temperatures, initializing the runs at the true valueof ␳˜ , with ␪ fixed as 10. The CPU times (on an 800-MHz Pentium
III processor) for data set 8 were ⵑ30 hr for infs and Recombine
and 30 sec for ␳ˆPAC-B.

I. Include all values of ␳ for which loge(LPAC-B(␳)) is
within 2 of the maximum.
II. Include all values of loge(␳) within ⫾1.96␴ of ␳ˆ PAC-B,
where ␴ is the square root of the inverse of minus
the second derivative (found numerically) of thelog of the PAC-B likelihood curve [as a function ofloge(␳)] evaluated at ␳ˆ PAC-A ⫽ ␳ˆ PAC-B.
The rationale for looking at such C.I.’s is that, understandard asymptotic theory for likelihood estimation,
C.I.’s constructed in this way using the true likelihood
curve would include the true value of ␳ ⵑ95% of thetime. (For I this follows from the asymptotic ␹2 distribution of the log-likelihood-ratio statistic; for II it follows
from asymptotic normal distribution of the MLE.)
Figure 6 shows the coverage properties for C.I.’s produced using the two methods (i.e., the proportion oftimes that C.I.’s formed using each method contained
the true value of ␳) for the data sets used to obtain
Figure 5c. For moderate sequence length both methods
produce C.I.’s that are slightly anticonservative, withcoverage properties that approach ⵑ0.91, comparedwith the expectation of ⵑ0.95 under asymptotic theory.
On the basis of these results we speculate that the curvature of the PAC-B-likelihood curve does not deviate
grossly from that of the true-likelihood curve. We notethat the coverage properties are also closer to asymptotic
expectations than those reported by Fearnhead andDonnelly (2002) for their composite likelihood using
the same methods of C.I. construction.
Comparison with other methods: We compared the LPAC-Blikelihood curves with likelihood curves obtained fromthree other methods: the full-data coalescent methodof FD (implemented in the computer program infs) andthe pairwise composite-likelihood methods of Hudson[2001; implemented by one of us (N.L.), using tables
available from R. Hudson’s website] and McVean et al.
(2002) (implemented in the computer program LDhat).

Linkage Disequilibrium and Recombination

2221

shown from infs are profile likelihoods for ␳ at the true
value of ␪; Hudson’s method and our method avoid
explicitly estimating ␪; LDhat estimates ␪ using an analog of Watterson’s estimate, but allowing for multiple
mutations.
Notwithstanding these issues, we attempt to draw
some general conclusions:Figure 6.—Empirical coverage properties of confidence
intervals produced using two different methods described in
the text. Each number is based on analysis of 1000 data sets
and shows the proportion of cases in which the C.I. contained
the true value of ␳ used to generate the data. The data sets
used are the same as those used to produce Figure 5c.

Figure 7 shows likelihood curves obtained using eachmethod for the 20 data sets considered by Wall (2000),
which were simulated under the standard coalescentmodel with ␳ ⫽ ␪ ⫽ 3.0 across a region of physicallength 1 and were kindly supplied by J. D. Wall. (These
likelihood curves are plotted with ␳ on the x-axis, rather
than log(␳), because infs and LDhat output likelihood
curves for evenly spaced values of ␳.)
Interpreting the results of this comparison is slightly
tricky. Unlike the other three methods we consider, thefull-data coalescent method can, in principle, provide
a fully accurate representation of the true-likelihood
curve. As such it is tempting to treat this as a “goldstandard” against which to compare the other methods.
However, as mentioned previously, even for the rather
small value of ␳ ⫽ 3 used to generate these data, accurate
approximation of the true-likelihood curve may be computationally impractical. Indeed, the estimated effective
sample sizes (ESSs) obtained for these data sets, shownat the top of each part of Figure 7, suggest that weshould not place much confidence in the accuracy ofmany of the curves. Our attempts to obtain more accurate likelihood curves by performing longer runs forsome of the data sets (numbers 15 and 16) actually produced smaller estimated ESSs, suggesting that the effective
sample sizes quoted for the other data sets are optimistic(see FD for further discussion of this problem). A further complication in comparing the methods is thatboth our method and that of McVean et al. (2002) allow(implicitly and explicitly, respectively) for the possibilityof multiple mutations, and thus the likelihoods fromthese methods are in some sense not directly comparable with those from the other two methods. Finally, we
note that the methods deal in different ways with the
unknown mutation parameter ␪: the likelihood curves

1. In general, the likelihood curves produced by thefour methods seem to agree rather more closely thanmight have been expected. (Compare, for example,
the variability here with the variability observed for
different runs of a single method in FD.) However,the closeness of the agreement between the methods
differs appreciably across data sets. Data set 12 consists of only four sites, three of which are singletons,
and so the differences in the curves for this data
set seem not to be particularly interesting. We wereunable to discern a systematic reason for the larger
differences among methods observed in some of theother data sets (e.g., 16).
2. The two pairwise composite-likelihood methods tend
to produce likelihood curves that are slightly more
peaked than those of the other two methods. Thismight be expected since, as pointed out by McVean
et al. (2002), pairwise composite-likelihood curves
are typically more peaked than the true-likelihood
curve because they treat each pair of sites as independent, when in fact many pairs are highly dependent.
3. The method implemented in LDhat, which allows
for multiple mutations, tends to achieve its maximum
at larger values of ␳ than does Hudson’s method,which does not allow for multiple mutations. This is
surprising; indeed, the opposite might have been
expected, since multiple mutations could be usedin place of recombination events to explain certain
patterns of LD. One possible explanation is that the
run lengths we used for computing the likelihoodin LDhat might be insufficient (we used the default
values).
4. Different orderings of the haplotypes can give PAClikelihood curves that differ appreciably from one
another. In addition, the maximum of the likelihood
curve based on the average over several orderings
tends to be toward the left end of the distributionof maxima obtained from different orderings. Thisis because, although not shown in Figure 7, the curves
with maxima at smaller values of ␳ tend to be larger(in absolute value) than those with maxima at larger
values of ␳ (presumably because they correspondto orderings of the haplotypes that, in some sense,
require fewer recombination events to explain them)and thus contribute more to the average. Although
this dependence on ordering is bothersome, in simulation studies (results not shown) we have found thatthe variability in the position of the maxima of the
PAC likelihood over different orderings of the haplo-

2222

N. Li and M. StephensFigure 7.—Comparison of the relative PAC-likelihood curves, with coalescent-based and pairwise composite relative-likelihood
curves, for the first 20 data sets in Wall (2000). In each case the relative likelihood is obtained by normalizing each likelihood
curve to have a maximum of 1. The light gray lines show 20 PAC-likelihood curves, each from a different random order of the
haplotypes, and the solid black line is based on the PAC likelihood averaged over the 20 random orders. The other lines
correspond to likelihood curves computed using the methods of: FD, implemented in the computer program infs (red dashedline); McVean et al. (2002), implemented in LDhat (blue dotted line); and Hudson (2001), using the table generated by program
eh written by Hudson (cyan dot-dashed line). The effective sample size (ESS) for infs at the MLE is given for each data set above
the graph and is a measure of the confidence infs has in its estimated likelihood curve (the larger the better). Results for infs
for all data sets except 15 and 16 were kindly provided to us by P. Fearnhead and were obtained using between 50,000 and
5,000,000 iterations. Results for data sets 15 and 16 were obtained by ourselves using 10,000,000 iterations.

types is typically small compared with the uncertainty
in estimation of ␳.
VARIABLE RECOMBINATION RATE

Models for variation in recombination rate: One ofour main motivations for developing this model is to
explore fine-scale variation in recombination rates. A
simple (no interference) model for variation in recombination rates is that crossovers in a single meiosis occur
as an inhomogeneous Poisson process, of rate c(x) at
position x. Here we consider two specific cases of thisgeneral model:
1. A simple single-hotspot model, wherec(x) ⫽␭c

冦c

for a ⱕ x ⱕ b,otherwise.

(4)

Here c represents the background rate of crossover,a and b represent the left and right ends of thehotspot region, and ␭ (⬎1) quantifies the magnitudeof the recombination hotspot. The PAC likelihoodfor this model is a function of four parameters: a, b,
␭, and ␳ ⫽ 4N c.
2. A more general model, where if x is a position between markers j and j ⫹ 1, thenc(x) ⫽ ␭j c.

(5)Here c represents the background rate of crossover,and ␭j is a multiplier controlling how the crossover
rate between markers j and j ⫹ 1 deviates from thebackground rate. The PAC likelihood for this model
is a function of the parameters ␭1, . . . ␭S⫺1 (where S
is the number of SNPs) and ␳ ⫽ 4 Nc.

Linkage Disequilibrium and Recombination

2223

TABLE 3
Performance of the simple single-hotspot model, for data sets simulated under various demographic scenarios,with a hotspot of magnitude ␭ ⫽ 10 and with no hotspot (i.e., ␭ ⫽ 1)
␭ ⫽ 10All sites
Common sites (f ⬎ 0.1)
Island (mixed)
Island (single)
Expansion (t ⫽ 500)
Expansion (t ⫽ 5000)

␭⫽1

Power

Med ␭ˆ

Med |Err|

Type I error

Med ␭ˆ

Med |Err|

0.90
0.81
0.92
0.94
0.94
0.53

13.02
16.57
11.35
10.24
24.29
11.66

0.19
0.29
0.18
0.16
0.40
0.36

0.04
0.05
0.07
0.07
0.13
0.07

0.95
1.13
0.90
1.01
1.41
1.23

0.45
0.57
0.33
0.31
0.56
0.70

In each case, the first column shows the proportion of “significant” LR tests (LLR ⬎ 1.92) for testing thenull hypothesis of no hotspot, the second column shows the median estimate of ␭, and the third column shows
the median of |Err(␭ˆ , ␭)|. Each number is based on results for 200 simulated data sets.

For the simple single-hotspot model it is straightforward to obtain numerically the maximum PAC-likelihood estimates for all four parameters simultaneously,
although in the examples that we consider we assume
that a and b are known and maximize the PAC likelihoodin terms of ␭ and ␳. The evidence for the presence ofa hotspot can be summarized by the log-likelihood ratio(LLR) for the null hypothesis of no hotspot, H0: ␭ ⫽ 1
vs. the alternative H1: ␭ ⬎ 1. If ␳0 denotes the value of␳ that maximizes LPAC-B under H0, and ␳1 and ␭1 denote the values of ␳ and ␭ that maximize LPAC-B underH1, then
LLR ⫽ loge LPAC-B(␳1, ␭1)/LPAC-B(␳0 , ␭ ⫽ 1),(6)

and large values of LLR represent evidence for the existence of a hotspot. Under standard asymptotic theory,
two times LLR would have (asymptotically) a chi-square
distribution on 1 d.f., and so rejecting H0 if LLR ⬎ 1.92
would give a hypothesis test with a type I error rate of
0.05. Although it seemed unlikely that standard asymptotic theory would apply here, we found that for data
sets simulated under the null hypothesis, rejecting H0
for LLR ⬎ 1.92 gave empirical type I error rates close
to 0.05 (Table 3), which provides some guidance as towhat might be considered a “large” value of LLR.
For the second, more general model, obtaining maximum PAC-likelihood estimates for the parameters creates problems. First, the maximum-likelihood estimates
are not unique (indeed there are infinitely many ofthem), because multiplying all the ␭j by any constant,and dividing ␳ by the same constant, gives exactly thesame likelihood. (In technical terms, the parameters
are said to be unidentifiable.) Second, even if the identifiability problem is solved (for example, by first obtaining an estimate for ␳ assuming that there is no hotspot and then fixing this when estimating the otherparameters) there is the practical problem that the likelihood curve for some ␭j will often be very flat, resulting
in estimates for many ␭j being very close to (or equalto) either 0 or infinity. This seems undesirable: if thelikelihood for a particular ␭j is very flat, this indicatesthat there is little information about the recombination
rate in that marker interval, in which case it seems sensible to estimate that the recombination rate is close to
the background rate (i.e., ␭j ⬇ 1), rather than (closeto) infinitely bigger or smaller!
To solve both these problems, we assume a “prior”
distribution for the ␭j’s: specifically that the ␭j’s are independent and identically distributed, with log10(␭j) ⵑ N(0,0.52). This prior was chosen to allow occasional deviations from the background rate of recombination by afactor of 10 or more (with probability ⵑ95%, ␭j liesin the range 0.1–10). This choice of prior could be
motivated from a Bayesian viewpoint as reflecting ourprior beliefs about the ␭j’s, but it also has the more
pragmatic justification that identifying variations of thiskind of magnitude seems both interesting and, perhaps,
attainable. We consider alternative prior specificationsin the discussion.
In principle, given the prior distribution for the ␭j’s
described above, we could also place a prior distribution
on ␳ and obtain an approximation to the posterior distribution of all parameters, using Markov chain Monte
Carlo, for example. Although this would be our preferred approach, for simplicity we avoid this here andinstead use the following ad hoc two-stage approach:first, obtain point estimates for ␳ and ␭ by jointly maximizing the product of LPAC(␳, ␭) and the prior densityof ␭; second, obtain a “posterior distribution” for each␭j by fixing all other parameters at their estimated values, discretizing the prior on ␭j (truncated at ␭j ⫽ ⫾3),and computing the corresponding discretized posterior
distribution as being proportional to the prior times
the PAC likelihood. For data sets with a large numberof sites the first stage (optimization over ␳, ␭) can bevery time consuming, requiring large numbers of evaluations of the likelihood function. Further, it seems unlikely that the simple optimization method we used will
reliably find the global maximum of the likelihood surface. Both these problems could be alleviated by ex-

2224

N. Li and M. Stephens

ploiting the fact that the derivatives of the PAC likelihood
can be computed efficiently, but we do not pursue thishere.
Power to detect recombination hotspots and robustness: In this section we assume that there is a single
recombination hotspot (Equation 4), whose putative
position is known, and examine the power of our model
to detect the hotspot under various assumptions about
the population demography and SNP marker ascertainment. Although the assumptions made here (in particular, that there is a single hotspot with known putative
position) are unrealistic, they provide a convenient
framework within which to examine quantitatively thepower of our approach and how it is affected by population demographic history and marker ascertainment
schemes.
We consider the following scenarios:
1. Constant-size randomly mating population, allmarkers
2. Constant-size randomly mating population, only
markers at frequency ⬎0.1
3. Exponentially expanding population, with expansion starting t ⫽ 500 generations ago
4. Exponentially expanding population, with expansion starting t ⫽ 5000 generations ago
5. Haplotypes sampled from a structured population,
consisting of two islands exchanging migrants at a
rate of one per generation (scaled migration parameter 4Nm ⫽ 4)
6. Haplotypes sampled from only one of the islands inthe structured population described above.
These last four models are the same as those considered in Pritchard and Przeworski (2001). In the two
expanding-population scenarios, the population was assumed constant sized until t generations ago, when it
started to expand exponentially, continuing until thepresent. The current population size N0 is set to be 105
and the population growth rate is chosen, as a functionof t, to match the expected diversity in a population of
constant size 104. (The necessary growth rates of ␣ ⫽1960, 350 for t ⫽ 500,5000 were kindly provided by M.
Przeworski.)For each scenario we simulated data sets under thesimple single-hotspot model described above, using
mksample and the postprocessing algorithm describedin Appendix C. For the first two scenarios each data
set was simulated to have ⵑ50 segregating sites and 60
chromosomes, with a ⫽ 0.4, b ⫽ 0.5, ␳ ⫽ 20, and ␭ ⫽10 [these values were chosen to approximately matchvalues for the TAP2 data from Jeffreys et al. (2000)
considered in Figure 9]. For the expanding population
scenarios we set ␳ ⫽ 4N0c ⫽ 200, and for the structured
scenarios ␳ ⫽ 20 within each population. For each scenario we also simulated data sets under the same conditions, but with no hotspot (i.e., ␭ ⫽ 1).
We applied the likelihood-ratio test to each data set

to test the null hypothesis H0: ␭ ⫽ 1 against the alternative H1: ␭ ⬎ 1 (Table 3). For the scenarios not involving
population expansion, the test gave type I error ratesof ⵑ0.05 when applied to data without a hotspot and apower of ⵑ0.90 when applied to data simulated with ahotspot, although the test based on just the common
SNPs had a slightly reduced power. The two scenarios
involving population expansion gave either a substantial
reduction in power or an inflated type I error rate
(which is in some sense equivalent to a reduction in
power). This might be due to a reduction in the numberof “common” SNPs under these scenarios, as common
SNPs tend to be most informative for estimating recombination rates.
We also examined the robustness of estimates of ␭under the various scenarios (Figure 3). As noted by FD,
the recombination rate ␳ ⫽ 4Nc depends on how theeffective population size N is defined. In contrast, thedefinition of the parameter ␭ does not depend on how
the effective population size is defined, and so we mighthope that estimation of ␭ will be robust to departuresfrom the assumption of a constant-sized panmictic population. For the levels of population structure we usedin our simulations this does indeed appear to be thecase—in both cases estimates were more accurate thanthose for the sample from a single random-mating population, perhaps because population structure makes recombinants easier to “spot.” As might be expected, estimates based only on common SNPs were less accurate
than those based on all SNPs. A drop in accuracy is also
evident for the scenarios simulated under population
expansion, probably again due to a reduction in thenumber of “common” SNPs under these scenarios.
Some of the scenarios also resulted in an upward bias for
estimates of ␭, notably one of the expansion scenarios inwhich the median of the estimates was almost 2.5 times
the true value.
Estimating recombination rates along a region: Simulated data: We fitted the more general varying recombination rate model to the simulated data used to produce
the pairwise LD plots in Figure 1; the results are shownin Figure 8. From the latter figure we might conclude,
correctly, that the data sets corresponding to the top left
and bottom middle parts had recombination hotspots
somewhere in the region 0.4–0.6. We might also conclude that the other data sets had no hotspots, which
would be correct except in the case of the bottom left
figure, which was actually generated from data with ahotspot between 0.4 and 0.5. One possible reason thatthe hotspot shows up less well in this case is that there
are fewer sites at high frequency (⬎0.15) in this data
set. Despite the fact that we might have been misled in
one case out of six, we view Figure 8 as considerably
more informative than Figure 1, from which we find it
difficult to draw any conclusions.
TAP2 data: Jeffreys et al. (2000) used patterns ofLD (measured by haplotype diversity) in a population

Linkage Disequilibrium and Recombination

2225

Figure 8.—Estimates of variation from background recombination rate within each marker interval for the same simulated
data sets that were used to produce Figure 1. Top left, bottom left, and bottom middle correspond to data sets simulated witha single hotspot of magnitude ␭ ⫽ 10 between positions 0.4 and 0.5. The other parts correspond to data simulated with constantrecombination rate across the region.

sample to refine the location of a putative recombination hotspot in the human TAP2 gene and provided amore detailed characterization of its properties throughsperm typing. The population sample consists of 30
individuals from the United Kingdom typed at 47 polymorphisms (45 SNPs, two insertion/deletions) across 9.7
kb, with haplotypes determined by allele-specific PCR.
Through analysis of sperm crossover events Jeffreys
et al. (2000) identified a region of increased crossover
intensity, located approximately in the interval from 4
to 5.2 kb.
We fitted the simple single-hotspot model to the haplotype data (kindly provided in convenient electronic
format by A. Jeffreys), assuming a hotspot between 4
and 5.2 kb, and obtained estimates of ␳ ⫽ 1.3/kb and␭ ⫽ 12 (95% C.I. [6,21]), with a LLR of 12, indicatingstrong evidence for the presence of the hotspot. Ourestimate for the average magnitude of the hotspot, ␭ ⫽12 times the background rate, agrees well with thesperm-typing results from Jeffreys et al. (2000). In particular, Jeffreys et al. (2000; their Figure 4) observed
128 crossovers within the interval 4–5.2 kb in 2.4 million
progenitor molecules, giving an average rate of 4.4 cM/
Mb, which is 11 times the approximate background rate(for males) that of 0.4 cM/Mb they quote, (although

they warn that this estimate of the background rateshould be “treated with caution”). Since our estimate
is based on a population sample, it is actually an estimateof the magnitude of the hotspot in the sex-averaged
crossover rates. Jeffreys et al. (2000) point out that thecrossover rate in this region appears to be substantiallyhigher in females than in males, and so the sex-averaged
crossover rates are likely to be dominated by the female
crossover process. Our results therefore suggest thatthe crossover rate within the 1.2-kb hotspot in female
meioses is roughly an order of magnitude higher thanthe (female) background rate.
Jeffreys et al. (2000) found that if one assumes thatthe sex-averaged background recombination rate isequal to the male background recombination rate of⬇0.4 cM/Mb (a number that we again emphasize theysuggest should be treated with caution), the observed
patterns of LD in the population sample appear consistent with an effective population size of N ⫽ 100,000,
which contrasts with the more commonly quoted number for humans of N ⫽ 10,000. Our estimate of ␳ ⫽1.3/kb supports their analysis, as it corresponds to N ⬇
84,000 if the background sex-averaged recombination
rate is assumed to be 0.4 cM/Mb. As pointed out byJeffreys et al. (2000), one possible explanation for this2226

N. Li and M. StephensFigure 9.—TAP2 data: MLE with posterior
quantiles. Results of fitting the more generalmodel for varying recombination rate to the TAP2
data from Jeffreys et al. (2000). The plot showsthe estimated value and 25, 10, and 5% posterior
quantiles for ␭j in each marker interval along thechromosome. The vertical lines indicate the approximate boundaries of the hotspot identified
by Jeffreys et al. (2000).

is differences between male and female recombination
rates—in particular, a sex-averaged background rate
across the 9.7-kb region of 3.4 cM/Mb would give N ⬇
10,000. An alternative (or additional) explanation, suggested to us by M. Przeworski (personal communication), is that gene conversion events not detected by
the sperm-typing experiments could partially accountfor the unusually large estimated effective population
size.
Figure 9 shows the estimates of ␭j and posterior quantiles obtained by fitting the more general model forrecombination rate variation to the TAP2 haplotype
data. The hotspot in the region 4–5.2 kb is fairly clear,with some suggestion that it may extend slightly farther
to the right than 5.2 kb. The peak of the recombination
hotspot is estimated as being ⵑ14 times the backgroundrate. In the interval corresponding to this peak the posterior probability that ␭j ⬎ 7 is 75% (compared with aprior probability of ⬍5%). However, the large numberof parameters estimated within this more general model
results in generally poor precision for the estimate ofeach ␭j. In particular, confidence in the estimates isprobably not sufficient to conclude that the three subpeaks present in Figure 9, within the hotspot, correspond to actual variation in the hotspot intensity.
Lipoprotein lipase data: The lipoprotein lipase (LPL)
data (Clark et al. 1998; Nickerson et al. 1998) consistof 9.7 kb of genomic DNA sequence from the human
lipoprotein lipase gene from 71 individuals from Jackson, Mississippi (n ⫽ 24), Rochester, Minnesota (n ⫽23), and North Karelia, Finland (n ⫽ 24). In the published data, the haplotypic phase for 69 sites was either
determined by experiment or estimated by Clark’s(1990) algorithm. Although the use of a statisticalmethod to infer some of the phases means that there
is some possibility that not all the published haplotypes
are completely correct, the majority seem likely to beaccurate, and in this analysis we assume them to be
known without error. On the basis of patterns of LD
and on the results of phylogenetic-based methods that
attempt to infer ancestral recombination events, Templeton et al. (2000) suggested the existence of a putative
recombination hotspot between [2987,4872].
Table 4 shows the results of fitting the simple singlehotspot model to the whole data set and to the datafrom each subpopulation individually, assuming a hotspot from 3 to 5 kb. Figure 10 shows the results forfitting the more general model for recombination ratevariation. Overall, these results seem to support theexistence of the putative hotspot, although there is considerable variation in the strength of the evidence (asmeasured by the LLR), and of the estimated magnitudeof the hotspot, in different subpopulations. We note thatthe apparent magnitude of the hotspot in the Finnish
population is smaller in Figure 10 than in Table 4, due
to the affect of the prior. There is also tentative evidence,

TABLE 4
Results of fitting the simple single-hotspot model to the
LPL data, to each subpopulation sample individually,
and to the combined sample
Population

␳/kb

␭ˆ

C.I.

LLR

JacksonRochesterFinland
Combined

7.5
0.74
0.15
7.0

2.7
12
104
6.5

[1.3, 4.8]
[5.0, 25]
[55, 183][4.1, 8.5]

3
8
22
23

The first two columns show estimated values for ␳ and ␭
assuming a hotspot between 3 and 5 kb along the sequence.
C.I. gives the range of values of ␭ whose log-likelihood is within
1.92 of the log-likelihood at ␭ˆ (for ␳ fixed at its estimatedvalue). LLR gives the value of the log-likelihood ratio for
testing the null hypothesis of ␭ ⫽ 1.

Linkage Disequilibrium and Recombination

2227

Figure 10.—Results of fitting
the more general model for varying recombination rate to the LPL
haplotype data from Clark et al.
(1998). The plot shows the estimated value and 25, 10, and 5%
posterior quantiles for ␭j in eachmarker interval along the chromosome. The vertical dashed lines
indicate the boundaries of the putative recombination hotspot identified by Templeton et al. (2000).

mostly from the Jackson sample, for a smaller-magnitude hotspot between 8 and 9 kb. Although no intervalin that region produces a very large estimate for ␭i, theclustering together of three intervals with moderate ␭iprovides stronger evidence than any one of these estimates taken separately.
Our results are consistent with those from Fearnhead
and Donnelly (2002), who found evidence for the[2987,4872] hotspot in the samples from Rochester andFinland, but not in those from Jackson. In addition,
both we and Fearnhead and Donnelly (2002) foundthat the Rochester and Finland samples give muchsmaller estimates for ␳ than the Jackson sample gives,probably reflecting smaller effective population sizes asa result of a recent bottleneck. One general advantageof the approach we take here, over the method of considering segments of the chromosome separately as do
Fearnhead and Donnelly (2002), is that it uses thepatterns of LD not only between markers within the
hotspot, but also between markers either side of the hotspot, to estimate the magnitude of the hotspot. This
may explain why we detected a signal (albeit a modestone) for a [2987,4872] hotspot in the Jackson sample,
where Fearnhead and Donnelly (2002) did not.
The large differences among estimates for ␭ from thethree separate population samples are surprising. To
examine whether this might be simply due to poor precision for these estimates in one or more of the populations (due, for example, to too much or too little diversity), we constructed ⵑ95% confidence intervals for ␭
using an analog of method I in the Construction of confidence intervals section (column C.I. in Table 4). Although

the coverage probabilities for these C.I.’s are unlikely
to be 95%, they give some indication of the curvatureof the likelihood surface, and the fact that not one of
the three intervals overlaps with either of the othertwo suggests that the hotspot intensity may indeed vary
among the three populations. Additional simulation results (not shown) suggest that the larger effective population size of the Jackson population should actually
increase power to detect the hotspot compared with the
Finnish population, and so differences in effective population sizes do not appear to explain our results. An
alternative explanation is the bias we observed for ourestimates of ␭ under certain expansion scenarios (Figure 3), which might partially explain the large estimateof ␭ in the Finnish population, for example, although
it is unclear whether this is enough to account for thefact that the estimated ␭ is almost 40 times greater thanthat in the Jackson population. Biological mechanismsthat could lead to different patterns of recombination
rate heterogeneity in different populations are known
to exist (e.g., Jeffreys and Neumann 2002), and thekinds of methods we introduce here should be helpful
in determining how often this occurs in practice.
DISCUSSION

In this article we have introduced a new statisticalmodel relating patterns of LD at multiple loci to theunderlying recombination rate and examined its effectiveness for inferring the underlying rate of recombination. Another potential application of our model is tomethods for LD (association) mapping in “case-control”

2228

N. Li and M. Stephensstudies, where chromosomes have been collected andtyped for both case and control individuals. Several authors, including McPeek and Strahs (1999), Morris
et al. (2000), and Liu et al. (2001), have developed methods to use genetic types at multiple loci to perform
association mapping for case-control studies. Thesemethods aim to improve on other common methods—
which typically test small groups of markers, one groupat a time, for association with a trait—by considering
data at many SNP markers simultaneously. Although
the methods differ in details, broadly speaking they all
pursue a strategy of assuming that (subsets of) the casechromosomes share some region identical by descent
about a causal mutation and as a result will be moresimilar than would be expected by chance. The challenge then is to identify regions where (subsets of) thecase chromosomes are more similar than would be expected by chance. Models of LD play a key role here,
because what would be expected “by chance” depends
critically on the amount of LD among loci. In particular,
correlations among loci will cause chromosomes to tend
to be more similar by chance than if the loci were independent. McPeek and Strahs (1999) use a first-orderMarkov chain to model LD, so that the probability ofobserving types (x1, . . . , xL) at L loci along a chromosome
is Pr(x1) Pr(x2|x1) Pr(x3|x2) . . . Pr(xL|xL⫺1), where theconditional probabilities Pr(xr|xr⫺1) are estimated using
the control chromosomes. This model was also adopted
by Morris et al. (2000). While the first-order Markov
assumption is better than assuming that the loci are
independent and may suffice if there is little LD amongmarkers, it seems not to be a good model for LD ingeneral. In particular, it fails to capture the fact thatmarkers may be in weak LD with neighboring markers,
but in strong LD with more distant markers. Although
McPeek and Strahs (1999) mention that higher-orderMarkov models might better model LD, such models
seem unlikely to be helpful in practice because of the
difficulty of estimating all the necessary parameters. The
model we have introduced here provides a parsimonious method for modeling LD: even the more generalmodel for varying recombination rates has fewer parameters than the first-order Markov model used previously.
Further, in these kinds of applications, where estimationof underlying recombination rates may be of only indirect interest, the usefulness of our model will depend
only on whether Pr(h1, . . . , hn|␳) is a sensible distribution
for h1, . . . , hn for some value of the parameters ␳, even if
this ␳ does not correspond precisely to the backgroundrecombination rate scaled by the effective population
size. Under these circumstances our two approximations ␲A and ␲B should perform almost identically, andso ␲A might be preferred on the grounds that it is simpler to understand and implement and is more amenable to theoretical study.
Another model for LD across multiple sites, introduced by Daly et al. (2001), is based on the empiricalobservation that in some regions of the genome LD
exhibits a “block-like” structure. Daly et al. (2001)
model each observed haplotype as a mosaic of “ancestral
haplotypes,” with the transition rates among these ancestral states (representing the “historical recombination
frequency” between each pair of consecutive markers)
being estimated by maximum likelihood. The ancestral
haplotypes are identified by an initial scan for regionsof low haplotype diversity, although in principle they
could instead be treated as parameters in the model.
Daly et al. (2001) used this model to produce a summary of patterns of LD that illustrates the haplotype
structure in their data more clearly, and in more detail,than would plots of pairwise LD measures. However, it
is currently unclear to what extent this model might behelpful for applications involving statistical inference,or prediction, particularly in regions where patterns ofLD are less block-like.
Several challenges might arise in applying our method
to real data, which we have ignored here. In particular,
we have assumed in our examples that haplotypes are
known and that there are no missing genotypes or genotyping errors. A new version of the software package
PHASE (Stephens et al. 2001) is under development,
and it will deal with these problems by incorporating
the PAC likelihood into a Markov chain Monte Carloalgorithm to jointly estimate the recombination rate
parameters, haplotypes, missing genotypes, and potential locations of genotyping errors. This algorithm also
produces a method for estimating haplotypes that takes
account of decay of LD along chromosomes. Preliminary results for simulated data suggest that these ideas
result in slightly more accurate haplotype estimates thandoes the method described in Stephens et al. (2001).
There are also biological aspects of real data that wehave not accounted for here, including, for example, geneconversion, whose effect on patterns of LD in humans
has been the subject of considerable recent interest (seeFrisse et al. 2001, for example). The effect that thepresence of gene conversion will have on our method
will vary, depending on how the tract length—about
which little is known in humans—compares with themarker density. Gene conversion events with very small
tract lengths compared to the marker density will only
rarely involve a typed marker and so will tend to havea small effect on our method unless such events are
extremely common. Conversely, gene conversion eventswith longer tract lengths—comparable to the typicaldistance between markers—will often affect one ormore markers and will tend to look like double-crossover events to our method. The presence of gene conversion with this kind of tract length will thus elevate ourestimates of recombination rate, perhaps substantially,
and regions with elevated rates of such gene conversion
may appear as recombination hotspots in our method.
In principle the PAC model could be extended to account explicitly for gene conversion by suitable modifi-

Linkage Disequilibrium and Recombination

cation of the conditional distribution ␲. A concretesuggestion for how to achieve this would be to augment
the space of the hidden Markov model for the mosaic
process (described in detail in appendix a) to include
both the current and previous “copied” chromosome
and then to modify the Markov jump process to makejumps back to the previously copied chromosome more
likely than jumps to other chromosomes. However, thiswould greatly increase the computational expense of
the model, making it unappealing in practice. A more
attractive possibility would be to settle for modeling only
those gene conversion events that affect a single marker(which, depending on tract length and marker density,
may be the vast majority of gene conversion events affecting patterns of LD). This would require only a simple
modification of the conditional distribution (it could be
handled similarly to the way that mutations are currently
handled), with essentially no increase in the computation required.
Another aspect of real data that we have not accounted for explicitly is population structure. Our simulation results in Figure 3 suggest that for the purposes ofidentifying recombination hotspots our method is robust
to a certain amount of population structure. Nevertheless,modeling population structure explicitly might prove
helpful in some settings. For example, it could be used
to extend methods for detecting population structurefrom unlinked markers (e.g., Pritchard et al. 2000) toallow them to be applied to sets of tightly linked markers. Again, a natural approach is to modify the conditional distribution ␲ to account explicitly for population
structure. One suggestion is to modify the copying process in the k ⫹ 1st chromosome (see appendix a) sothat, rather than being equally likely to copy all r existing
chromosomes, it is more likely to copy chromosomesfrom the same population than chromosomes from adifferent population. This would effectively model population structure by increasing the probability of seeing
similar chromosomes in the same population, comparedwith in different populations. We are currently investigating the effectiveness of a similar idea for LD mappingin case-control studies: treating cases and controls as separate populations and examining whether there appears
to be evidence in some regions for the case chromosomes
to be more similar to other case chromosomes than tocontrol chromosomes.
While we have concentrated here on models for biallelic loci, the ideas we have introduced could also be
used to model LD among multiallelic loci such as microsatellites. There is a natural analog of ␲A for loci with
K alleles (see also the conditional distribution for
K-allele loci suggested in FD), and this could form astarting point for further investigation.
To deal with the problem that the PAC likelihood
depends on the order in which the haplotypes are considered, we have chosen to average the likelihood over
several random orders. One possible alternative would

2229

have been to use the pseudo-likelihood (Besag 1974)
based on our conditional distribution,n

L pseudo(␳) ⫽ 兿 ␲(hk |H⫺k),

(7)

k⫽1

where H⫺k denotes the set of all haplotypes excluding
hk. The pseudo-likelihood, by definition, does not depend on the ordering of the haplotypes. This idea ismore along the lines of the way that these conditional
distributions are used in Stephens et al. (2001). However, in preliminary studies we found that this pseudolikelihood performed poorly for estimating ␳. Our intuitive explanation for this is that the pseudo-likelihoodin effect contains only information on the recombination that is occurring in the tips of the trees and noton the structure of the tree as a whole. (Interestingly,
under our approximation the first two haplotypes contain no information on ␳, so in some sense the information on ␳ comes from intermediate haplotypes.) Nonetheless, it is possible that the pseudo-likelihood may
prove useful in settings where estimating ␳ is not of
direct interest.
We have introduced here two models for variation inrecombination rate: a simple single-hotspot model anda more general model that allows recombination rates
to vary along the chromosome. Each of these models
has weaknesses. The simple single-hotspot model makessome unrealistic assumptions: the background recombination is unlikely to be constant and neither is therecombination rate within the hotspot; furthermore,
there could be more than one hotspot. The more general model makes few assumptions and allows more
flexible investigation of patterns of recombination ratevariation along a region. However, this extra flexibility
comes at the expense of the introduction of extra parameters, which can result in a reduction in the precision with which parameters can be estimated. Whenusing the model as a general model for LD, rather thanfor parameter estimation as we have concentrated onhere, the precision of parameter estimates may be unimportant, and the few assumptions made by the moregeneral model make it particularly attractive in this situation. When estimation of recombination rates is themain goal, the more general model may be viewed as mostsuited to exploratory data analysis, identifying plausible
positions for hotspots, whose magnitudes might be estimated by a more parsimonious model. In this situation
it might prove fruitful to consider modifying the moregeneral model by putting a more informative prior on
the ␭j’s. In particular, a prior in which the ␭j’s are correlated along the chromosome (e.g., an autoregressive
prior) would reduce the variance of parameter estimates atthe expense of assuming that changes in recombination
rate occur more or less smoothly along the chromosome(which may or may not be the case).
In assessing our model as a method for estimating
recombination rates from sequence data over moderate

2230

N. Li and M. Stephensgenomic regions, perhaps the most natural comparisons
to make are with the composite-likelihood methods of
Hudson (2001) and Fearnhead and Donnelly (2002).
(While some other methods based on summaries of thedata might be competitive with these approaches when
the recombination rate is assumed constant, they seem
likely to suffer from loss of information when fitting
models with more parameters, such as either of ourmodels for recombination rate variation.) Of the twocomposite-likelihood methods, although both are tractable for estimating ␳ over large genomic regions, only
Hudson’s method is comparable with our own in termsof computational expense: our method and Hudson’s
method typically take seconds, or less, to compute alikelihood, while Fearnhead and Donnelly’s method cantake hours per likelihood computation. Although thetime and effort expended in collecting these kinds ofdata make it not unreasonable to wait hours or days forthe results of analysis, the extra computational burden
may make Fearnhead and Donnelly’s method difficult
to extend to more general settings involving missing
genotype data, genotyping error, and/or unknown
haplotypic phase, for example. The approach of splitting the sequence data into contiguous segments also
has the disadvantage, noted earlier, of estimating recombination rates in a region only on the basis of sites
within the region and not sites either side of the region,
resulting in potential loss of information. Our limited
comparisons with Hudson’s method suggest that it performs similarly to our method for estimating the recombination rate when it is assumed constant across theregion. In principle, Hudson’s method could also be
applied to fit models of varying recombination rate
along the sequence, and the existence of more than
one method to fit such models would be welcome. Bothapproaches seem to offer considerable advantages overother available methods for modeling LD and inferring
patterns of recombination rate heterogeneity.

SOFTWARE

The methods discussed in this article have been implemented in a C⫹⫹ package, Hotspotter, which is freely
available for download at http://www.stat.washington.edu/
stephens/software.html.
We thank Peter Donnelly, Paul Fearnhead, Gil McVean, Jonathan
Pritchard, and Molly Przeworski for helpful conversations and Paul
Fearnhead, Richard Hudson, and Mary Kuhner for helpful advice
and/or results from software. Two anonymous referees gave helpfulcomments on the submitted manuscript. This work was supported bya grant from the University of Washington and National Institutes of
Health grant no. 1RO1HG/LM02585-01.

LITERATURE CITED
Besag, J., 1974 Spatial interaction and the statistical analysis of latticesystems. J. R. Stat. Soc. Ser. B 36: 192–236.

Clark, A. G., 1990 Inference of haplotypes from PCR-amplified
samples of diploid populations. Mol. Biol. Evol. 7 (2): 111–122.
Clark, A. G., K. M. Weiss, D. A. Nickerson, S. L. Taylor, A.
Buchanan et al., 1998 Haplotype structure and population genetics inferences from nucleotide-sequence variation in human
lipoprotein lipase. Am. J. Hum. Genet. 63: 595–612.
Daly, M. J., J. D. Rioux, S. F. Schaffner, T. J. Hudson and E. S.
Lander, 2001 High-resolution haplotype structure in the human genome. Nat. Genet. 29: 229–232.
Ewens, W. J., 1972 The sampling theory of selectively neutral alleles.
Theor. Popul. Biol. 3: 87–112.
Fearnhead, P. N., and P. Donnelly, 2001 Estimating recombination rates from population genetic data. Genetics 159: 1299–1318.
Fearnhead, P., and P. Donnelly, 2002 Approximate likelihood
methods for estimating local recombination rates. J. R. Stat. Soc.
Ser. B 64: 657–680.
Frisse, L., R. R. Hudson, A. Bartoszewicz, J. D. Wall, J. Donfack
et al., 2001 Gene conversion and different population histories
may explain the contrast between polymorphism and linkage
disequilibrium levels. Am. J. Hum. Genet. 69: 831–843.
Griffiths, R. C., and P. Marjoram, 1996 Ancestral inference fromsamples of DNA sequences with recombination. J. Comput. Biol.
3: 479–502.
Hammer, M. C., T. Karafet, A. Rasanayagam, E. T. Wood, T. K.
Altheide et al., 1998 Out of Africa and back again: nested
cladistic analysis of human Y chromosome variation. Mol. Biol.
Evol. 15 (4): 427–441.
Harding, R. M., S. M. Fullerton, R. C. Griffiths, J. Bond, M. J. Cox
et al., 1997 Archaic African and Asian lineages in the genetic
ancestry of modern humans. Am. J. Hum. Genet. 60: 772–789.
Hey, J., and J. Wakeley, 1997 A coalescent estimator of the population recombination rate. Genetics 145: 833–846.
Hudson, R. R., 1983 Properties of a neutral allele model with intragenic recombination. Theor. Popul. Biol. 23: 183–201.
Hudson, R. R., 1987 Estimating the recombination parameter of afinite population model without selection. Genet. Res. 50: 245–
250.
Hudson, R. R., 2001 Two-locus sampling distribution and their application. Genetics 159: 1805–1817.
Hudson, R. R., 2002 Generating samples under a Wright-Fisher
neutral model of genetic variation. Bioinformatics 18: 337–338.
Hudson, R. R., and N. L. Kaplan, 1985 Statistical properties of thenumber of recombination events in the history of a sample of
DNA sequences. Genetics 111: 147–164.
Jeffreys, A. J., and R. Neumann, 2002 Reciprocal crossover asymmetry and meiotic drive in a human recombination hot spot. Nat.
Genet. 31: 267–271.
Jeffreys, A. J., A. Ritchie and R. Neumann, 2000 High resolution
analysis of haplotype diversity and meiotic crossover in the human
TAP2 recombination hotspot. Hum. Mol. Genet. 9: 725–733.
Jeffreys, A. J., L. Kauppi and R. Neumann, 2001 Intensely punctate
meiotic recombination in the class II region of the major histocompatibility complex. Nat. Genet. 29: 217–222.
Johnson, G. C. L., L. Esposito, B. J. B. A. N. Smith, J. Heward, G. D.
Genova et al., 2001 Haplotype tagging for the identification ofcommon disease genes. Nat. Genet. 29: 233–235.
Kingman, J. F. C., 1982 The coalescent. Stoch. Proc. Appl. 13: 235–
248.
Kruglyak, L., 1999 Prospects for whole-genome linkage disequilibrium mapping of common disease genes. Nat. Genet. 22: 139–144.
Kuhner, M. K., J. Yamato and J. Felsenstein, 2000 Maximum likelihood estimation of recombination rates from population data.
Genetics 156: 1393–1401.
Liu, J. S., C. Sabatti, J. Teng, B. J. Keats and N. Risch, 2001 Bayesian analysis of haplotypes for linkage disequilibrium mapping.
Genome Res. 11: 1716–1724.
McPeek, M. S., and A. Strahs, 1999 Assessment of linkage disequilibrium by the decay of haplotype sharing, with application tofine-scale genetic mapping. Am. J. Hum. Genet. 65: 858–875.
McVean, G., P. Awadalla and P. Fearnhead, 2002 A coalescentbased method for detecting and estimating recombination fromgene sequences. Genetics 160: 1231–1241.
Morris, A. P., J. C. Whittaker and D. J. Balding, 2000 Bayesianfine-scale mapping of disease loci, by hidden Markov models.
Am. J. Hum. Genet. 67: 155–169.
Nickerson, D. A., S. L. Taylor, K. M. Weiss, A. G. Clark, R. G.

Linkage Disequilibrium and Recombination
Hutchinson et al., 1998 DNA sequence diversity in a 9.7-kb
region of the human lipoprotein lipase gene. Nat. Genet. 19:
233–240.
Nielsen, R., 2000 Estimation of population parameters and recombination rates from single nucleotide polymorphisms. Genetics
154: 931–942.
Olivier, M., V. I. Bustos, M. R. Levy, G. A. Smick, I. Moreno etal., 2001 Complex high-resolution linkage disequilibrium andhaplotype patterns of single-nucleotide polymorphisms in 2.5 Mbof sequence on human chromosome 21. Genomics 78: 64–72.
Press, W. H., S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1992 Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, New York.
Pritchard, J. K., and M. Przeworski, 2001 Linkage disequilibriumin humans: models and data. Am. J. Hum. Genet. 69: 1–14.
Pritchard, J. K., M. Stephens and P. Donnelly, 2000 Inferenceof population structure using multilocus genotype data. Genetics
155: 945–959.
Rabiner, L. R., 1989 A tutorial on hidden Markov models and selected applications in speech recognition. Proc. IEEE 77: 257–
286.
Stephens, M., and P. Donnelly, 2000 Inference in molecular population genetics. J. R. Stat. Soc. Ser. B 62: 605–655.
Stephens, M., N. J. Smith and P. Donnelly, 2001 A new statisticalmethod for haplotype reconstruction from population data. Am.
J. Hum. Genet. 68: 978–989.
Templeton, A. R., A. G. Clark, K. M. Weiss, D. A. Nickerson, E.
Boerwinkle et al., 2000 Recombination and mutation at the
LPL locus. Am. J. Hum. Genet. 66: 69–83.
Ueberhuber, C. W., 1997 Numeric Computation: Methods, Software,and Analysis, Vol. I. Springer-Verlag, Berlin-Heidelberg, Germany.
Wakeley, J., 1997 Using the variance of pairwise differences to estimate the recombination rate. Genet. Res. 69: 45–48.
Wall, J. D., 2000 A comparison of estimators of the population
recombination rate. Mol. Biol. Evol. 17: 156–163.
Wang, N., J. M. Akey, K. Zhang, R. Chakraborty and L. Jin, 2002
Distribution of recombination crossovers and the origin of haplotype blocks: the interplay of population history, recombination,and mutation. Am. J. Hum. Genet. 71: 1227–1234.
Communicating editor: S. Tavare

APPENDIX A:
THE CONDITIONAL DISTRIBUTION ␲ A

Here we give a formal description of ␲A. We also provide some additional motivation for the form of thisapproximation and describe briefly some of the variations on this form with which we have also experimented.
Formal description of ␲A: Let h1, . . . , hn denote then sampled haplotypes typed at S biallelic loci (SNPs).
Typically h1, . . . , hn would come from a sample of n
haploid individuals or n/2 diploid individuals. We assume that the distribution of the first haplotype is independent of ␳ [e.g., all 2S possible haplotypes are equally
likely, so ␲A(h1) ⫽ 1/2S]. Consider now the conditional
distribution of hk⫹1, given h1, . . . , hk, for k ⱖ 1. Recall(Figure 2) that hk⫹1 is an imperfect mosaic of h1, . . . ,hk. That is, for k ⱖ 1, at each SNP, hk⫹1 is a (possibly
imperfect) copy of one of h1, . . . , hk at that position.
Let Xj denote which haplotype hk⫹1 copies at site j (soXj 僆 {1,2, . . . , k}). For example, for haplotype h4A in
Figure 2, (X 1,X 2,X 3,X 4,X 5) ⫽ (3,3,2,2,2). To mimic theeffects of recombination, we model the X j as a Markov
chain on {1, . . . , k}, with Pr(X1 ⫽ x) ⫽ 1/k (x 僆 {1,. . . , k}), and

2231

Pr(X j⫹1 ⫽ x⬘|X j ⫽ x)⫽exp(⫺␳jd j/k) ⫹ (1 ⫺ exp(⫺␳ jd j/k))(1/k)

冦(1 ⫺ exp(⫺␳ d /k))(1/k)

if x⬘ ⫽ x;otherwise ,

j j

(A1)
where d j is the physical distance between markers j andj ⫹ 1 (assumed known); and ␳j ⫽ 4Nc j , where N is theeffective (diploid) population size, and cj is the averagerate of crossover per unit physical distance, per meiosis,
between sites j and j ⫹ 1 (so that cj d j is the genetic
distance between sites j and j ⫹ 1). This transition matrix
captures the idea that, if sites j and j ⫹ 1 are a small
genetic distance apart (i.e., cjd j is small), then they arehighly likely to copy the same chromosome (i.e., X j⫹1 ⫽X j). We note the following special cases used in thisarticle:1. Constant recombination rate (see estimating constant recombination rate): cj ⫽ c for all j.
2. Simple single-hotspot model (see Power to detect recombination hotspots and robustness): c j ⫽ c if markers j
and j ⫹ 1 are both outside the hotspot, and c j ⫽ ␭c
if markers j and j ⫹ 1 are both inside the hotspot.
(For brevity we omit details of the more tedious,
though straightforward, case in which one marker isin the hotspot and the other outside the hotspot.)
3. General variable recombination rate model (see Estimating recombination rates along a region): cj ⫽ ␭j c.
To mimic the effects of mutation, the copying process
may be imperfect: with probability k/(k ⫹ ␪˜ ) the copy
is exact, while with probability ␪˜ /(k ⫹ ␪˜ ), a “mutation”
will be applied to the copied haplotype. Specifically, if
hi,j denotes the allele (0 or 1) at site j in haplotype i,then, given the copying process X1, . . . , X S, the alleleshk⫹1,1, hk⫹1,2 , . . . , hk⫹1,S are independent, withPr(hk⫹1,j ⫽ a|Xj ⫽ x,h1, . . . , hk)⫽k/(k ⫹ ␪˜ ) ⫹ (1/2) ⫻ ␪˜ /(k ⫹ ␪˜ ),冦(1/2) ⫻ ␪˜/(k ⫹ ␪˜),h x, j ⫽ ahx , j ⬆ a .

(A2)[The factor of (1/2) appears in both cases, so that as␪˜ → ∞ both alleles become equally likely.]
We fix the value of ␪˜ to be
␪˜ ⫽

n⫺1

⫺1

冢 兺 m1 冣 ,

(A3)

m⫽1

where n is the total number of sampled haplotypes. (See
Motivation and Variation below for more discussion.)
Computation: Computing ␲A(hk⫹1|h1, . . . , hk) requires
a sum over all possible values of the X j , which can bedone efficiently using the forward part of the forwardbackward algorithm for hidden Markov models (e.g.,
Rabiner 1989). Specifically, let hk⫹1,ⱕj denote thetypes of the first j sites of haplotype hk⫹1, and let ␣j(x) ⫽Pr(hk⫹1,ⱕ j ,X j ⫽ x). Then ␣1(x) can be computed directly

2232

N. Li and M. Stephensfor x ⫽ 1, . . . k, and ␣2(x), . . . . ␣S(x) can be computed
recursively usingk␣j⫹1(x) ⫽ ␥j⫹1(x) 兺 ␣j (x⬘)Pr(X j⫹1 ⫽ x|X j ⫽ x⬘)

(A4)

x⬘⫽1

冢

冣

1 k
兺 ␣j (x⬘) ,k x⬘⫽1⫽ ␥j⫹1(x) pj ␣j(x) ⫹ (1 ⫺ pj)(A5)

where ␥j⫹1(x) ⫽ Pr(hk⫹1,j⫹1|X j⫹1 ⫽ x, h1, . . . , hk) is givenin (A2), and pj ⫽ exp(⫺␳j d j/k).
The value of ␲A(hk⫹1|h1, . . . , hk) can then be computedusingk␲A(hk⫹1|h1, . . . , hk) ⫽ 兺 ␣S(x).

(A6)

x⫽1

The second term in the parentheses of (A5) does notdepend on x and needs to be computed only once
for each j (as noted in FD). Thus the computational
complexity of ␲A(hk⫹1|h1, . . . , hk) increases linearly inthe number of SNPs and linearly in k. As a result, thecomputation of LPAC-A is linear in the number of SNPs
and quadratic in the number of chromosomes in thesample.
Motivation and variations: Although it seems intuitively sensible that the transition matrix in Figure 8
should have the property that the rate at which jumps
occur in the copying process should increase with ␳ anddecrease with the number of previous sampled chromosomes k, it is perhaps not so obvious why we chosethe rate ␳/k. Indeed, the empirical results in Figure 4
suggest that this rate is not quite ideal, and appendix
b describes a corrected rate, based on these empiricalresults. However, we can get some idea of why ␳/n is asensible starting point for the rate parameter from thefollowing informal argument. Assume that h1, . . . , hk⫹1
are a random sample of haplotypes from a neutrally
evolving, randomly mating, constant-sized population,and consider the unknown genealogical tree relating
h1, . . . , hk⫹1 at a single site. It follows from the Ewens
sampling formula that in this tree, the probability thathk⫹1 is separated by at least one mutation from each of
h1, . . . , hk (unconditional on the actual values of h1,. . . , hk) is ␪/(k ⫹ ␪), where ␪ ⫽ 4N␮, and ␮ is theprobability of mutation per meiosis at that site. Similarly,
if we consider marking on the tree recombination events
that occur between this site and the next site, the probability that there will be at least one such event separatinghk⫹1 from each of h1, . . . , hk is ␳/(k ⫹ ␳), where ␳ ⫽4Nc and c is the probability of recombination between
two adjacent sites per meiosis. Since ␳ is small, ␳/(k ⫹␳) ⬇ ␳/k, giving the rate that we used. (We emphasizethat this is not intended to be a formal argument andin particular that it is unclear how our mosaic process
relates formally to the genealogical tree relating the
haplotypes. It is merely intended to provide additional
motivation for the use of this rate and perhaps to stimulate research into a more formal connection.)

The reason for our choice of ␪˜ is that ␪˜ 兺nm⫺⫽111/m is the
expected number of mutation events at a single site on
the genealogical tree relating a random sample of nchromosomes, so (A3) gives a priori an expected numberof mutation events at each site of 1 (although it does
not force the number of mutations to be exactly 1,and so our method should be somewhat robust to thepresence of multiple mutations at some sites).
We performed simulation experiments along the linesof those used to produce Figure 4 to see whether variations on the conditional distribution ␲A described above
might eliminate the bias we observed for ␲A. In particular, we tried using values for ␪˜ that were up to four times
bigger or smaller than that in (A3); estimating ␪˜ fromthe data; replacing the transition probability in (A1)with a transition probability of ␳di/(k ⫹ ␳d i), as in FD;
and making use of the more complex mutation mechanism (involving Gaussian quadrature) used in Stephensand Donnelly (2000) and in FD. Although these different variations gave different quantitative results, they
all produced similar qualitative patterns, and in particular the bias we observed for ␲A remained in every variation that we tried. We therefore resorted to the empirical correction described in appendix b below.

APPENDIX B:␲ B, A BIAS-CORRECTED VERSION OF ␲ A

To correct the bias observed in the results for ␳ˆ PAC-A,
we modified the transition matrix in (A1) by replacing
␳j by ␦j␳ j, where ␦j ⫽ exp(a ⫹ b log10␳j). The intercept
a and slope b are interpolated on the basis of the numberof haplotypes n and segregating sites S in the data (Figure 1), using tensor product interpolation with natural
cubic splines first in the direction of varying n and thenin the direction of varying S (Ueberhuber 1997).

APPENDIX C: SIMULATING DATA WITH A
RECOMBINATION HOTSPOTWe use the following algorithm to postprocess theoutput from mksample (Hudson 2002) to simulate data
under the simple single-hotspot model for recombination variation. Suppose we would like to simulate a sample with approximately S segregating sites. The background recombination rate is ␳. A hotspot of width w ⫽(b ⫺ a) lies between positions a and b, with recombination rate ␭␳ where ␭ ⬎ 1. We follow these steps:
1. Simulate samples with S⬘ ⫽ (1 ⫹ w(␭ ⫺ 1))S segregating sites and constant recombination rate ␳⬘ ⫽ (1 ⫹w(␭ ⫺ 1))␳.
2. Multiply the position of each site by a factor of 1 ⫹w(␭ ⫺ 1) so that the total length of the haplotypes

Linkage Disequilibrium and Recombinationis 1 ⫹ w(␭ ⫺ 1) instead of 1 (and the backgroundrecombination rate is ␳).
3. For sites within a and a ⫹ w␭, randomly delete themwith probability 1 ⫺ 1/␭.
4. For the remaining sites within a and a ⫹ w␭, shrink
the distance of adjacent sites by a factor of ␭.
5. Shift the positions of the sites to the right of the2233

hotspot [subtract w(␭ ⫺ 1)] so that the total lengthis again 1.
Shrinking the distance between sites in the hotspot
produces the effect of elevated recombination rate. Deleting some sites keeps the mutation rate constant over
the region.