Am. J. Hum. Genet. 76:449–462, 2005

Accounting for Decay of Linkage Disequilibrium in Haplotype Inference
and Missing-Data Imputation
Matthew Stephens and Paul Scheet
Department of Statistics, University of Washington, Seattle

Although many algorithms exist for estimating haplotypes from genotype data, none of them take full account ofboth the decay of linkage disequilibrium (LD) with distance and the order and spacing of genotyped markers. Here,
we describe an algorithm that does take these factors into account, using a flexible model for the decay of LD with
distance that can handle both “blocklike” and “nonblocklike” patterns of LD. We compare the accuracy of thisapproach with a range of other available algorithms in three ways: for reconstruction of randomly paired, molecularly determined male X chromosome haplotypes; for reconstruction of haplotypes obtained from trios in an
autosomal region; and for estimation of missing genotypes in 50 autosomal genes that have been completely
resequenced in 24 African Americans and 23 individuals of European descent. For the autosomal data sets, ournew approach clearly outperforms the best available methods, whereas its accuracy in inferring the X chromosome
haplotypes is only slightly superior. For estimation of missing genotypes, our method performed slightly better
when the two subsamples were combined than when they were analyzed separately, which illustrates its robustness
to population stratification. Our method is implemented in the software package PHASE (v2.1.1), available fromthe Stephens Lab Web site.

Introduction
At autosomal loci, the genetic material carried by a diploid individual can be thought of as being composed of
two haplotypes, each containing the genetic informationfrom one of the two homologous chromosomes. Knowledge of the haplotypes carried by sampled individuals
would be helpful in many settings, such as linkage-disequilibrium (LD) mapping or attempting to make inferences regarding evolutionary mechanisms, such as selection or recombination, that may be acting on a region.
However, although technologies continue to develop,
current high-throughput approaches to genotyping do
not provide haplotype information. Although, in somestudies, collection of data from related individuals may
allow haplotypes to be inferred, in general, such data
may be costly or impossible to collect. These factors havegenerated considerable interest in computational andstatistical approaches for inferring haplotypes from unphased genotype data in population samples.
All computational and statistical approaches to haplotype inference exploit LD, which is the nonrandom
association of alleles among linked loci. LD tends to
Received August 9, 2004; accepted for publication January 4, 2005;
electronically published January 31, 2005.
Address for correspondence and reprints: Dr. Matthew Stephens, University of Washington, Department of Statistics, Box 354322, Seattle,
WA 98195-4322. E-mail: stephens@stat.washington.edu
䉷 2005 by The American Society of Human Genetics. All rights reserved.
0002-9297/2005/7603-0009$15.00

decay with distance; that is, there tends to be less LD
between loci that are far apart than between loci thatare close together. Earlier approaches to haplotype estimation—which include the algorithm of Clark (1990),
maximum likelihood via the EM algorithm (Excoffierand Slatkin 1995; Hawley and Kidd 1995; Long et al.
1995), and the Bayesian approach of Stephens et al.
2001—ignore this fact, in that haplotype estimates fromthose approaches are independent of locus spacing andeven of locus order.
As far as we are aware, the first haplotype-inference
method to produce results that can depend on locus
order was introduced by Niu et al. (2002). Those authors introduced the idea of “partition ligation” (PL),which approaches data sets containing several loci byfirst dividing them into segments containing a small
number (∼8) of contiguous loci and then inferring haplotypes within each segment before iteratively combining results from adjacent segments. This framework was
subsequently adopted by others, including Qin et al.
(2002), Stephens and Donnelly (2003), and Lin et al.
(2004), for different methods for inferring haplotypes
within each segment and for combining segments. The
most obvious benefits of PL are the reduction of computing times and the increase in the size of data set thatcan be tackled effectively. However, a side effect is that
estimates of phase at any locus will depend more on
data at nearby loci than on data at more-distant loci.
Although this effect is presumably beneficial, it is unfortunately diluted by the fact that all existing approaches
449

450

ignore the order and spacing of loci within segments
formed.
Other recent methods for haplotype inference (e.g.,
Greenspan and Geiger [2004], Halperin and Eskin[2004], Kimmel and Shamir [2004]) also involve dividing data into segments of consecutive loci. However,unlike the methods described above, the division into
segments is made with reference to observed patternsof LD. Motivated by the observation that some regionsof the human genome appear to contain “blocks” of
low haplotype diversity, these methods attempt to ensure that the segments formed correspond, in some sense,
to “blocks” of high LD and of low haplotype diversity.
Excoffier et al. (2003) take a different approach, using
a window of neighboring loci when estimating phase at
any given position, with window widths allowed to vary
according to local levels of LD. None of these methods
take into account locus spacing or the decay of LD
within the blocks or windows formed.
Here, we describe a new algorithm for haplotype inference that modifies the approaches of Stephens et al.
(2001) and Stephens and Donnelly (2003), to take explicit account of locus spacing and the decay of LD with
distance. These previous approaches are Bayesian andmake use of an “approximate coalescent” prior to reflect the fact that haplotypes in a population tend togroup together in clusters of similar haplotypes. Our
new algorithm can be thought of as modifying this prior
to the “coalescent with recombination” (Hudson 1991),
by use of recently developed models for how patternsof LD among multiple loci depend on the underlying
recombination rate (Fearnhead and Donnelly 2001; Li
and Stephens 2003). The assumptions underlying thiscoalescent process include the fact that the population
has been evolving with a constant size for a long timeand is randomly mating (so is in Hardy-Weinberg equilibrium) and that the locus under study is evolving neutrally. Since none of these assumptions will hold in realstudies, it may be helpful to think of our prior not as
an attempt to perform inference under a formal model
but rather as an attempt to capture salient qualitative
features of real data—notably, the tendency for haplotypes to cluster together (because of shared ancestry)
and the fact that which haplotypes cluster together may
change as one moves along the sequence (because ofrecombination). Since the underlying recombination
rate is typically unknown and may vary on a fine scale(Crawford et al. 2004; McVean et al. 2004), we adopt
a flexible model that estimates a different recombination
rate in each SNP interval. In this way, we allow for both
gradual decay of LD with distance and for more-abrupt
decay of LD, as might be observed across a recombination hotspot, for example. As a result, the method isequally applicable, whether or not the region understudy exhibits “blocklike” patterns of LD.

Am. J. Hum. Genet. 76:449–462, 2005

We compare the accuracy of this new algorithm withseveral of the numerous algorithms for haplotype inference now available (table 1). First, we assess accuracyof haplotype estimates, using X chromosome haplotypes determined from males (Lin et al. 2002). Second,
we compare a subset of the methods on the haplotype
data, determined from genotypes in trios in an autosomal region by Daly et al. (2001). Last, we assess theaccuracy with which haplotype-inference methods areable to predict missing genotypes, using 50 genes resequenced by the University of Washington–Fred Hutchinson Cancer Research Center (UW-FHCRC) Variation
Discovery Resource (Carlson et al. 2003). A considerable advantage of this last approach to comparing methods is that it requires only unphased genotype data,
which are relatively abundant. Since real data almost
invariably contain some missing genotypes, the accuracy with which such genotypes can be predicted fromother (nonmissing) genotypes is of direct practical interest. Further, accuracy in predicting missing genotypes
provides another test of how accurately the assumptions
underlying each method capture patterns of real haplotype variation and, thus, an indirect way to assess
which methods are likely to provide the most accurate
haplotype estimates. Accurate methods for imputing
missing data may also be helpful in developing methods
for LD mapping and choosing haplotype-tagging SNPs
(Johnson et al. 2001) (see the “Discussion” section).
Despite these considerations, there appear to be few
previous published comparisons of the accuracy of different methods for this task.
All our assessments suggest that haplotype-inference
methods based on an approximate coalescent prior are
consistently more accurate than are the other methods
we consider. Of coalescent-based approaches, our new
approach clearly outperforms the previous approachesof Stephens et al. (2001) and Stephens and Donnelly(2003) (which ignore the decay of LD with distance) for
the autosomal data sets, whereas its performances ininferring the X chromosome haplotypes is only slightly
superior.
Methods
We now give an overview of the algorithm and discuss
associated theoretical issues. A more detailed description
is given in appendix A.
Our algorithm is similar to those used by Crawford
et al. (2004) and Ptak et al. (2004), although the focus
here is different: in those articles, the unknown underlying haplotypes were treated as “nuisance parameters,”
and the underlying recombination process was the object of interest, whereas here we treat the recombination
process as a “nuisance parameter” and assess the accuracy of estimated haplotypes. Another difference isTable 1
Summary of the Haplotype Inference Algorithms for Unrelated Individuals
Algorithm
PHASE v2 recom (-MR)

Description

Bayesian method with approximate “coalescent with recombination” prior, capturing the fact thateach sampled haplotype tends to be similar to another haplotype or to a mosaic of other
haplotypes.
PHASE v2 hybrid (-MQ)
Hybrid algorithm reduces computing time over -MR by assumption of no recombination for the majority of the computation, before incorporating recombination for the final steps.
PHASE v2 no recom (-MS) Bayesian method with approximate “coalescent without recombination” prior, on the basis of theidea that sampled haplotypes look similar to other haplotypes. Ignores decay of LD with distance.
Bayes-Dirichlet
Bayesian method with Dirichlet prior on population haplotype frequencies. Ignores similarity of haplotypes. (Our own implementation of the algorithm in Haplotyper, mentioned below.)
Arlequin 3.0a (ELB)
Bayesian method. Prior takes some account of similarity of haplotypes. Bases inference for each locus
on data in a window of nearby loci, which is allowed to vary in size on the basis of local levels ofLD.
Haplotyper
Bayesian method with a Dirichlet prior on the haplotype frequencies. Ignores similarity of haplotypes.
PL-EM
A maximum-likelihood approach. Uses EM algorithm plus computational trick (PL) to obtain estimates of haplotype frequencies and uses these estimates to find most probable haplotypes for eachindividual.
snphap v1.2
Similar to PL-EM above, but uses different tricks to reduce computational demands.
hap
Based on an ad hoc modification of the Dirichlet prior for population haplotype frequencies.
hap2
An improved implementation of hap above.
HAP
Algorithm partitions loci into blocks of low haplotype diversity; assumes that haplotypes withinblocks will conform approximately to a “perfect phylogeny” (i.e., no recombination or repeat
mutation).

Reference(s)Present study

Present study
Stephens et al. 2001; Stephens and Donnelly 2003
Stephens et al. 2001; Niu et al. 2002
Excoffier et al. 2003

Niu et al. 2002
Qin et al. 2002

David Clayton Software Web page
Lin et al. 2002
Lin et al. 2004
Eskin et al. 2003a, 2003b; Halperin and Eskin 2004

452

Am. J. Hum. Genet. 76:449–462, 2005

that, here, we adopt a more flexible model for the underlying recombination process.
Let G p (G1 , … ,Gn ) denote the observed genotypesof the n individuals at L loci; let H p (H1 , … ,Hn ) denote the actual (unobserved) haplotypes, where Hi is
the pair of haplotypes of individual i. Note that Gi may
be “missing” alleles at some loci, in which case H includes estimates of the unobserved alleles. Let r p(r1 , … ,rL⫺1) denote the (unknown) vector of recombination rates between each pair of consecutive loci,scaled by the effective population size. That is,rl p

4Ne cl,dl

where Ne is the (unknown) diploid effective population
size, cl is the (unknown) probability of recombination
per generation between markers l and l ⫹ 1, and dl isthe physical distance between markers l and l ⫹ 1, assumed known. Informally, rl dl measures the expected
breakdown in LD across the interval spanned by markers l and l ⫹ 1, with large values corresponding to alarge expected breakdown in LD.
Our aim is to estimate H from the observed genotype
data G, taking account of information in the data—
particularly in the rate of decay of LD with distance—on the underlying recombination rates r. We do this bydeveloping a Markov chain–Monte Carlo (MCMC) algorithm to sample from the conditional distribution of
the haplotypes and recombination parameters, given the
genotype data Pr (H,rFG). Using this algorithm, we canobtain a point estimate for H (e.g., the estimated posterior mode for the haplotype pair of each individual),
together with measures of confidence in the accuracy of
individual haplotype estimates.
In outline, our algorithm starts with initial guesses
for H and r and a random ordering of the individuals
n and then iterates the following steps many times:
1. for each individual i in turn (in the order given by
the ordering n), update the individual’s pair of haplotypes by sampling Hi from Pr (HiFGi ,H⫺i ,r), where H⫺iis the set of current guesses for the haplotypes of all
individuals except i;
2. propose a new value for r and accept it or rejectit according to the Metropolis-Hastings (MH) acceptance probability (eq. [2]); and3. propose a new value for n (by proposing to switch
two randomly chosen individuals) and accept it or rejectit according to the MH acceptance probability.

Step 1 can be thought of as estimating the haplotypes
for individual i, taking into account the individual’s
genotype data (Gi), the haplotypes of all other individuals (H⫺i), and the underlying patterns of LD (r). Several
previous haplotype-inference algorithms (e.g., those of

Stephens et al. [2001], Niu et al. [2002], and Stephensand Donnelly [2003]) adopt a similar structure, although
without the parameter r, and accomplish this step indifferent ways, depending on the underlying assumptions. Here, we make use of the conditional distributions from Fearnhead and Donnelly (2001) and Li andStephens (2003), which, in estimation of haplotypes inindividual i, favors haplotypes that are similar to a mosaic of the haplotypes possessed by other individuals.
Break points in such a mosaic can occur anywhere but
are most likely to occur in marker intervals that arelarge or have a higher recombination rate (higher rl);see figure 1 for illustration. Thus, in inferring the haplotypes, the method takes into account marker spacing
and the decay of LD with distance.
More precisely, we use
Pr (Hi p {h,h  }FGi ,H⫺i ,r ) ∝ (2 ⫺ dhh )# p(hFH⫺i ,r)p(h FH⫺i ,r) ,(1)

where dhh  p 1 if h p h  and p 0 otherwise; the conditional distribution p is a modification of that of Fearnhead and Donnelly (2001), as described below. (Thesecond term on the right side should be “p(h FH⫺i ,h,r),”but we omitted “h” from this conditioning for convenience of implementation and for the happy side effectof producing an expression that is symmetrical in h andh , which, for our choice of p, is not otherwise
guaranteed.)
Step 2 of the algorithm can be thought of as estimating recombination rates from patterns of LD in thecurrent estimated haplotypes. The MH-acceptance probability is given by[

A p min 1,

Q(r  r r)L(r  )p(r  ),Q(r r r  )L(r)p(r)](2)

where Q(r r r  ) is the probability of proposing a moveto r , given that we are currently at r, p(r) is a prior onr, and L(r) p Pr (HFr) is the likelihood for r. We use
the PAC-B likelihood of Li and Stephens (2003) for
L(r). (This likelihood depends on the order in which
the haplotypes are considered; we use the order given
by n to order the individuals and order the two haplotypes within each individual alphabetically.) Our prior
on r is similar to the “general recombination variation”
model used by Li and Stephens (2003). Specifically, we
assume that ri p r̄l i, where the prior on l i is thatthey are independent and identically distributed, withlog 10 l i normally distributed, with mean 0 and SD 0.5,
and the prior distribution for r̄ is uniform on a log scale,with the constraint 10⫺8 ! r̄ ! 10 3. Here, r̄ can bethought of as the background recombination rate acrossthe region and l i the factor by which ri deviates from

Stephens and Scheet: Accounting for Decay of Linkage Disequilibrium

453

Figure 1

Illustration of how p(hk⫹1Fh1, … ,hk) builds hk⫹1 as an imperfect mosaic of h1, … ,hk. The figure illustrates the case k p 3 andshows three possible values (h4A , h4B , and h4C ) for h4 , given h1 , h2 , and h3 . Each column of circles represents a SNP locus, with blackened andunblackened circles representing the two alleles. Each possible h4 can be thought of as having been created by “copying” parts of h1, h2 andh3. The shading in each case shows which haplotype was “copied” at each position along the chromosome and indicates whether the haplotype
is most closely related to h1 , h2 , or h3 . Changes in the shading along a haplotype represent ancestral recombination events; these are more likely
to occur between SNPs that are farther apart (e.g., SNPs 2 and 3 or 4 and 5) or that have a higher rate of recombination between them. The
imperfect nature of the copying process is exemplified at the third and fourth locus, where h4B has the blackened allele despite having “copied”
h2, which has the unblackened allele. The occurrence of several such “imperfections” on a single chunk indicate that the haplotype is relatively
highly diverged from the one that it copied in that region (see text).

this background. The prior on l i says that 95% of intervals between SNPs will have a recombination ratethat is within a factor of 10 of the background rate.
The range of allowed values for r̄ covers several ordersof magnitude on either side of the average value forhumans.
Step 3 of the algorithm is included partly to ensurethat the algorithm (provided it is run long enough) will
not depend on the individuals’ entry order in the input
file and partly to provide an order over which to compute the PAC-B likelihood in step 2. The approach issimilar to that taken by Stephens and Donnelly (2003),
except that they randomly ordered the individuals every
iteration, whereas we use an MH proposal (with acceptance probability computed using a uniform prioron all possible orderings and the PAC-B likelihood).
The above outline omits many details. Perhaps themost important of these is that strict implementation of
step 1 is not computationally tractable for data sets witha moderate number of loci, because it requires expression (1) to be evaluated for all haplotype pairs {h,h  }consistent with genotypes of individual i, and the number of such pairs may be huge. To avoid this problem,
we make use of PL (Niu et al. 2002; Stephens and Donnelly 2003), wherein the data are first divided into small
segments of consecutive loci, and the algorithm is thenapplied to each segment in turn before iteratively combining the results from adjacent segments. Each time
the algorithm is applied to a segment, a list is made of

the haplotypes that may “plausibly” occur (in at least
one individual) within that segment. This list is typicallyshort, compared with the total number of possible haplotypes. When two segments are combined, only haplotype pairs consistent with haplotypes on the corresponding lists are considered possible, and the fact thatthese lists are typically short keeps the computationalburden manageable. Other details of the algorithm aregiven in appendix A.
Our algorithm, in common with algorithms of Stephens et al. (2001) and Stephens and Donnelly (2003),
contains a “pseudo-Gibbs” step (Heckerman et al. 2000),in that the conditional distributions p used in our step
1 are defined directly, rather than being derived froman explicit prior distribution on H. Care must be takenwith the use of algorithms containing such steps, since
they are not, in general, guaranteed to converge to astationary distribution. In previous work, we have defended this approach and, in particular, have shown thatthe similar algorithms of Stephens et al. (2001) andStephens and Donnelly (2003) are guaranteed to converge to a stationary distribution because they are defined on a finite discrete space (see the appendix of Stephens and Donnelly [2003]). Unfortunately, this argument does not apply to this new algorithm, because
step 2 includes a continuous parameter, r. Theoreticalconvergence of this algorithm, therefore, remains an
open question; however, in our experience with a wide
range of examples—both real data and simulations—454

we have seen no indication of convergence problems
beyond those usually associated with these kinds of
MCMC schemes, and, in our comparisons, the algorithm produces, on average, more-accurate haplotype
estimates (by two different measures) than does the algorithm of Stephens and Donnelly (2003), which isguaranteed to converge.
Another theoretically dubious aspect of the algorithm
is our use of different conditional distributions (thoseof Fearnhead and Donnelly [2001] and Li and Stephens[2003]) in steps 1 and 2 of the algorithm; the reason
for this is explained below. We note that it would bepossible to use the models of Li and Stephens (2003)
to produce MCMC schemes that do not suffer fromeither of these theoretical weaknesses, but they would
require substantially more computational time and seem
to us unlikely to produce practically important gains in
performance.
Formal Description of p
The conditional distribution p used to compute expression (1) is a slight modification of the one introduced
by Fearnhead and Donnelly (2001) and is also similar
to the conditional distribution pA of Li and Stephens(2003). For simplicity, we change notation slightly, andlet h1 , … ,hk denote the (not necessarily distinct) haplotypes of k chromosomes sampled from a population, andlet hk⫹1 denote the haplotype of another chromosome
sampled from the same population. The conditional distribution p(hk⫹1Fh1 , … ,hk ,r) specifies the probability ofobserving a given value for hk⫹1, given the observed
h1 , … ,hk and the values for the underlying recombination parameters r. It is based on the idea that hk⫹1 will
look similar to a mosaic of h1 , … ,hk (fig. 1). Think ofhk⫹1 as being, at each locus l p 1, … ,L, a (possibly imperfect) copy of one of h1 , … ,hk, and let Xl denote which
haplotype hk⫹1 copies at locus l (so Xl 苸 {1,2, … ,k}). In
addition, we associate with the copying event at locus la time Tl, which affects the probability that hk⫹1 will beidentical to the allele it copied at locus l. Informally, one
can think of Xl as indicative of which of the haplotypes
h1 , … ,hk is most closely related to hk⫹1 at locus l andTl (or, more precisely, Tl /k) as indicative of how closely
related the two are. We allow Tl to take one of two possible values, which can be thought of as corresponding
to “closely” or “distantly” related. Specifically, Tl cantake values (t1 ,t 2 ) p (0.586,3.414), with respective marginal probabilities (w1 ,w2 ) p (0.854,0.146). These numbers were chosen as a very rough approximation of acontinuous exponential distribution with mean 1, on thebasis of Gaussian quadrature with Laguerre polynomials(Evans 1993); see also the appendix of Stephens andDonnelly (2000). It is in the introduction of these Tj
values that this conditional distribution differs fromAm. J. Hum. Genet. 76:449–462, 2005

those of Li and Stephens (2003) and is similar to theone from Fearnhead and Donnelly (2001).
To mimic the effects of recombination, we model{(Xl ,T):lp 1, … ,L} as a Markov chain on {1, … ,k}, withl
Pr (X1 p x,T1 p tr) p (1/k)wr for (x,r) 苸 {1, … ,k} #{1,2}, andPr (Xl⫹1 p x ,Tl⫹1 p trFXl p x,Tl p tr) p{

exp (⫺rl dl /k)⫹[1 ⫺ exp (⫺rl dl /k)](1/k)wr if x  p x and r  p r ;[1 ⫺ exp (⫺rl dl /k)](1/k)wr
otherwise .
(3)

The idea here is that jumps in the mosaic process occur
with probability 1 ⫺ exp (⫺rl dl /k) (and so are most likely
to occur in marker intervals in which rl dl is large) and,
when a jump occurs, it is equally likely to jump to copying any of the chromosomes. Note that the jump probability we use differs slightly from that of Fearnhead andDonnelly (2001), who use rl dl /(k ⫹ rl dl). The numericalvalues of the two expressions are similar for small (typical) values of rl dl, but our expression has the desirable
theoretical property that probabilities are unchanged by
the addition of a locus, with a missing allele in hk⫹1,
anywhere along the sequence.
To mimic the effects of mutation, we allow that theallele of hk⫹1 at locus l may not be an exact copy of
the allele that it “copied.” Specifically, we assume thatthe allele at locus l is the result of applying a number(m, which may be 0) of mutations to the allele that it
copied. We assume that m has a Poisson distribution
˜ /k, where v˜ is a parameter the value ofwith mean vTlwhich controls the frequency of mutations. (Thus, consistent with the interpretation of Tl as a measure ofdivergence, mutations are more likely at loci where Tl
is large.) The mutation mechanism at each locus l is
specified by a matrix Pl , assumed to be known, whose(i,j)th element is the probability that an offspring is oftype j, given that the progenitor is of type i and that amutation occurs. Thus, if hi,l denotes the allele at locus lin haplotype i, then, given the copying process (X,T) p(X1 , … ,XL ,T1 , … ,TL ), the alleles hk⫹1,1 ,hk⫹1,2 , … ,hk⫹1,L
are independent, withPr (hk⫹1,l p aFXl p x,Tl p t,h1 , … ,hk ,r)

冘
⬁

p

˜ m
(vt/k)m˜
exp (⫺vt/k)(P)hx,la .
m!
mp0

(4)

We approximate this infinite sum by the first 50 terms(which need be done only once for each locus, and theresults tabulated).
In all the examples we consider here, the data consistof biallelic SNP loci. For these, we use, as did Stephens

455

Stephens and Scheet: Accounting for Decay of Linkage Disequilibrium

et al. (2001), ṽ p ( 冘 mp1 m1 )⫺1 and a mutation process(P) in which a mutation at a locus always leads to achange from the current allele to the alternative allele.
In software, we have also implemented a generalized
stepwise mutation mechanism for microsatellite loci anda parent-independent mutation mechanism for other multiallelic loci, with ṽ estimated as was done by Stephenset al. (2001 [see their appendix A]).
Note that when rl p 0 for all l (i.e., no recombination), the conditional distribution p described above simplifies to the conditional distribution used by Stephenset al. (2001).
Computation of p requires a sum over all possible
values of (X,T):n⫺1

Pr (hk⫹1 p hFh1 , … ,hk ,r)
p

冘写
X,T

j

Pr (hk⫹1,l p aFXl ,Tl ,h1 , … ,hk ,r) Pr (X,T) ,(5)

which can be accomplished efficiently by use of standard
computational methods for hidden Markov models (asshown in Li and Stephens [2003]; e.g., their appendix
A). The computational effort required to compute p bythis approach increases linearly with the product of thenumber of loci and number of individuals.
Introducing T doubles the amount of computational
effort required to compute each conditional distribution,compared with the distribution used by Li and Stephens(2003), but we believe it to improve phasing accuracy.
For example, suppose that, at two closely linked SNPloci, we observe many homozygotes (e.g., genotypes AA
and AA) and one double heterozygote (e.g., genotypes
AC and AC). We would argue that the most likely explanation is that one of this individual’s two haplotypes
is relatively highly diverged from the others in the sample
and that the singleton mutations therefore lie together
on this haplotype; that is, that this last individual most
likely has one AA haplotype and one CC haplotype.
Indeed, this is the solution favored by most haplotypeinference methods, including maximum likelihood via
the EM algorithm, Clark’s algorithm, and previous versions of PHASE. However, without the use of the times
T, our algorithm would randomize the phase of thesesingletons, independently at each site, and therefore consider the alternative solution (haplotypes AC and CA)
equally plausible. With the Ts, our algorithm will identify one of the haplotypes as being highly diverged fromthe others and therefore tend to put the singletons all
together on the same haplotype, as we argue it should.
To what extent this improves the accuracy of phasing
for nonsingleton SNPs is unclear.

Results
Our new algorithm is implemented with the software
package PHASE v2 (v2.0 and subsequent versions include variants on the algorithm we describe here; resultsin this article are based on v2.1.1). The algorithm comesin two “flavors,” which we label “recom” and “hybrid,”or “-MR” and “-MQ,” after the switches used in thesoftware to invoke them. The difference between themis that the first uses the “coalescent with recombination”
at all times, whereas the second is a hybrid between thisnew recom algorithm and the algorithm of Stephens andDonnelly (2003) (in that it uses the coalescent without
recombination until the final merge in the PL, when it
switches to using the coalescent with recombination).
The hybrid algorithm can be considerably faster and, as
we demonstrate here, sacrifices relatively little in accuracy of haplotype estimates—although we caution that
we have not assessed accuracy of recombination-rate
parameter estimates for the hybrid algorithm. We compare the accuracy of these algorithms with the coalescent-based method of Stephens and Donnelly (2003),
which does not take explicit account of the decay of LD
and which we term “no recom” (invoked with the -MS
switch in PHASE v2), and with several other available
algorithms for haplotype inference (table 1).
Haplotype Inference: X Chromosome Data
Our first comparison uses SNP data analyzed and described by Lin et al. (2002). The data consist of X chromosome haplotypes derived from 40 unrelated males.
The haplotypes comprise eight regions, which have arange of 87–327 kb and include 45–165 segregatingsites. For each of the eight genes, we created 100 data
sets, each consisting of 20 pseudoindividuals, created byrandomly pairing the 40 chromosomes (as in previous
comparisons using these data [Lin et al. 2002; Stephensand Donnelly 2003]).
We ran all computer programs except Arlequin 3.0awith their default or recommended settings, except as
noted below. (Results for Arlequin 3.0a were kindly
supplied by L. Excoffier.) For Haplotyper and PL-EM,“rounds” was set to 20. HAP was accessed via the HAP
Webserver (accessed October 2003 and December 2003).
For snphap, we reduced the lower posterior-trimming
threshold to 0.001, in an attempt to obtain haplotype
reconstructions for all individuals in the sample. The
input format for hap2 (kindly provided by S. Lin) allows
for the inclusion of family-level data; we input the individuals as unrelated parents with no offspring, and weused the following additional settings: 10,000 iterations,5,000 burn-in, 20 thinning, (0, 0.5) minor-allele frequency limits, and FD F blocks with 0.8 threshold.
To evaluate the accuracy of the haplotype reconstruc-456

Am. J. Hum. Genet. 76:449–462, 2005

tions for the various methods, we calculated two error
rates:
1. the individual error rate, defined as the proportionof ambiguous individuals whose haplotypes are not completely correct, ignoring in each individual any positions
with missing data, and2. the switch error, which measures the proportion ofheterozygote positions whose phase is incorrectly inferred relative to the previous heterozygote position.
This is 1 minus the switch accuracy defined by Lin etal. (2002).
For haplotypes with large numbers of SNPs or that coverlarger genetic distances, it is difficult to correctly infer
the entire haplotypes, so, in these cases, the switch error seems a more informative criterion for comparing
methods.
Results for all methods are given in table 2. Results
for the method of Stephens and Donnelly (2003) differ
slightly from those given in that article, partly because

we used a different 100 random pairings of the haplotypes and partly because we corrected errors in the locus
order used by Stephens and Donnelly (2003) for two of
the genes (GLA and TRP). For each gene, the rankingof the algorithms is generally similar for both the error
measures we use. For six of eight genes, one of our new
coalescent-based methods has the lowest individual-error rate; for seven of eight genes, one of them has thelowest switch-error rate. Next best is the coalescentbased method without recombination, followed by the
Bayes-Dirichlet method (which is our implementation ofan algorithm similar to the one underlying Haplotyper)and the method of Excoffier et al. (2003) implemented
in Arlequin 3.0. One possible explanation for the superior performance of Bayes-Dirichlet compared with
Haplotyper on these data is that the version of Haplotyper we used appeared to have problems accurately
imputing missing data, which may have adversely affected its accuracy in some cases (although we excluded
positions with missing data when scoring the algorithms).

Table 2Comparison of Accuracy of Methods for Reconstructing Haplotypes
INDIVIDUAL-ERROR RATE FOR GENE
ALGORITHM
PHASE v2 recom (-MR)
PHASE v2 hybrid (-MQ)
PHASE v2 no recom (-MS)
Bayes-Dirichlet
Arlequin (ELB)
hap2
hap
PL-EM
snphap
Haplotyper
HAP

GLRA2

MAOA

KCND1

ATR

GLA

TRPC5

BRS3

MECP2

Average (Ⳳ SE)

.73
.76
.77
.82
.73
.78
.79
.82
.82
.89
.93

.52
.53
.53
.58
.56
.62
.61
.63
.60
.76
NA

.46
.47
.48
.56
.53
.72
.54
.70
.68
.72
.86

.47
.47
.48
.57
.57
.58
.62
NA
.74
.72
.89

.70
.69
.70
.76
.76
.88
.89
NA
.90
.79
.98

.66
.63
.63
.61
.74
.79
.58
NA
NA
.72
NA

.65
.65
.66
.74
.62
.70
.72
.76
.74
.79
.82

.77
.78
.79
.84
.77
.84
.85
NA
NA
.64
.92

.62 (Ⳳ.015).62 (.015).63 (.015).68 (.015).66 (.013).74 (.013).70 (.016).73a (.014).75a (.017).75 (.009).90a (.009)

SWITCH-ERROR RATE FOR GENE

PHASE v2 recom (-MR)
PHASE v2 hybrid (-MQ)
PHASE v2 no recom (-MS)
Bayes-Dirichlet
Arlequin (ELB)
hap2
hap
PL-EM
snphap
Haplotyper
HAP

GLRA2

MAOA

KCND1

ATR

GLA

TRPC5

BRS3

MECP2

Average (Ⳳ SE)

.08
.09
.10
.11
.11
.09
.14
.12
.16
.16
.18

.07
.06
.06
.08
.09
.07
.10
.09
.10
.12
NA

.14
.14
.14
.17
.22
.20
.22
.23
.24
.27
.35

.17
.17
.18
.25
.23
.23
.29
NA
.36
.32
.41

.11
.11
.12
.14
.14
.14
.22
NA
.23
.16
.31

.17
.15
.15
.14
.22
.19
.13
NA
NA
.20
NA

.12
.10
.10
.14
.20
.11
.14
.14
.14
.15
.18

.16
.16
.17
.20
.20
.17
.23
NA
NA
.19
.31

.12 (.005).12 (.005).13 (.005).15 (.007).18 (.007).15 (.007).18 (.008).15a (.015).21a (.015).20 (.008).29a (.015)

NOTE.—The best performance in each column is highlighted in bold italics. The eight genes are the data sets on which methods
were tested by Lin et al. (2002). The individual- and switch-error rates are defined in the text. Estimated SDs for the average error
rate (final column) are given in parentheses (although, for assessment of the significance of an observed difference between two
methods, it would be best to take account of the paired nature of the data by computing an SE for the difference in error rates).
The data for Haplotyper and hap are taken from Lin et al. (2002). For two of the genes, indicated by “NA,” and for unknown
reasons, HAP would not process the input files we submitted. Similarly, for four genes, also indicated by “NA,” PL-EM andsnphap failed to complete runs for the majority of simulations.
aAverage computed with partial data.

457

Stephens and Scheet: Accounting for Decay of Linkage Disequilibrium

Another possible factor is that our implementation of
PL used a different criterion for deciding how long thelist of “plausible” haplotypes should be. Our criterion
may tend to lead to longer lists, which should increase
accuracy at the expense of extra computation. Overall,
hap and hap2 exhibit performances slightly superior tothose of PL-EM and Haplotyper, which are comparablewith each other.
One unfortunate aspect of data sets formed from randomly pairing X chromosomes is that they can include
genotypes in which one allele—and only one allele—isknown to be missing. This does not typically occur inreal genotype data, and it may be thought to adversely
affect relative performance of some of the methods, particularly HAP and snphap, since the versions of these
programs we used were unable to use information on
the single observed allele at such loci. To investigate thispossibility, we reran PHASE-MR, ignoring the single observed allele at such loci. The results were almost identical to those obtained when those alleles were not ignored, suggesting that this unusual pattern of missing
alleles is probably not the explanation for the relatively
poor performance of HAP for these data. (It is worth
noting that the comparisons of Halperin and Eskin[2004], in which HAP was competitive with PHASE,
used an earlier version of PHASE, implementing the algorithm of Stephens et al. [2001]).
Haplotype Inference: Autosomal Data
A slightly surprising feature of our results is that our
new algorithms appear to produce only small gains in
performance over coalescent-based algorithms that do
not explicitly model decay in LD. Since LD decays more
quickly on the autosomes than on the X chromosome(Schaffner 2004), it is possible that relative performanceof the methods might differ on autosomal data. We
therefore applied our coalescent-based algorithms to the
genotype data for the children in the 129 trios describedby Daly et al. (2001) and compared estimated haplotypes with the actual haplotypes inferred from the parental genotypes (ignoring genotypes for which phase
could not be so inferred) (kindly supplied to us in aconvenient electronic form by G. Kimmel). Since Halperin and Eskin (2004) also used these data for comparisons, we included HAP in this comparison, using theresults posted on the HAP Webserver.
For these data, our new algorithms (PHASE-MR and-MQ) produced more-accurate haplotype estimates thaneither the coalescent-based algorithm that ignored recombination (PHASE-MS) and HAP: respective individual error rates were 0.42, 0.39, 0.50, and 0.51; switch
error rates were 0.030, 0.034, 0.043, and 0.057. Clearly,
the relative performance here of HAP is better than for
the X chromosome data: the switch-error rate is less thandouble that of our new approach, compared with morethan double (on average) in the X chromosome comparisons. This might have been expected, since patterns
observed in these data partly motivated some of the assumptions underlying HAP.
Missing Data
To assess the accuracy of the algorithms for imputing
missing alleles, we used genotype data of 24 African
Americans (AA) and 23 individuals of European descent(ED) at 50 genes that we randomly selected from theUW-FHCRC Variation Discovery Resource Web site
(September 2003). These genes have been completely
resequenced and contain 15–230 segregating sites, withan average of 85 per gene. For convenience, we ignored
the few triallelic SNPs and multisite insertion-deletion
polymorphisms in the data.
To assess how well the various haplotyping methods
impute missing genotype data, we introduced into thedata 5% artificial missingness (in addition to the 4.6%
native missingness). We used two different patterns of
missingness: missing alleles and missing genotypes. In
the first, each observed allele was denoted as missing
with probability 0.05, independent of all other genotypic
data and missingness patterns. In the second, each observed genotype was denoted as missing with probability
0.05, again independently. Neither of these patterns islikely to capture all aspects of patterns of missingness
in real data, in which variations in DNA quality or molecular effects can cause some individuals and some sites
to have more than their fair share of missing data. For
this reason, absolute accuracy of methods of this test
may tend to be better than for real data. However, we
would expect relative performances to be similar.
To provide a baseline against which to compare methods, we implemented a naive “straw man” approach to
imputing the genotypes. For the alleles-missing pattern,
this method imputes the most common allele in the sample at that site; for the genotypes-missing pattern, it imputes the most common genotype in the sample at thatsite. Thus, the straw man ignores LD when imputing
missing data, and the size of the improvement of eachmethod over the straw man indicates the effectivenesswith which that method exploits patterns of LD in its
imputation algorithm.
To measure the accuracy of imputed genotypes, weused the genotype-imputation error rate, which we define as the total number of imputed genotypes not identical to the original genotype, divided by the total number of imputed genotypes (both totals computed acrossall 50 genes, ignoring native missing genotypes).
We analyzed the 50 genes with each method, using
identical missingness patterns for each method, in two
different ways. First, we analyzed all 47 individuals to-458

Am. J. Hum. Genet. 76:449–462, 2005

gether. Second, we analyzed the AA and ED samples
separately.
Table 3 shows the accuracy of the genotypes estimatedby each method. Although Arlequin (ELB) can be used
to analyze data sets with missing data, it does not estimate the missing alleles and so was omitted from thiscomparison. For the majority of genes, PL-EM and Haplotyper failed to present resolved haplotypes for at least
one of the three samples (AA, ED, and combined sample), which made it difficult to summarize results for
those methods. We omitted results for Haplotyper but
included the results for PL-EM for the eight genes thatit successfully processed in all three samples. Similarly,snphap successfully processed only 34 of the 50 genes.
HAP, hap2, and snphap appeared not to deal fully withgenotypes for which only one allele is missing, and so
we omitted them from the “alleles missing independently” set of comparisons.
Our new algorithms outperformed all other methods,
including the PHASE no recom method. Further examination revealed that this improvement is achieved bya consistent improvement across many genes rather thana large improvement in a small number of genes. For
example, when analyzing the “genotype missing” data,with both samples together, the Bayes-Dirichlet method
produced more-accurate results than did our new method
for only 4/50 genes, and HAP produced more-accurate
results for 1/50 genes.
For all of the PHASE methods (and some of theothers), the accuracy of imputed genotypes was better
when analyzing the ED and AA samples together thanwhen analyzing them separately, although differenceswere small in absolute terms. On the basis of this, we
speculate that, for these data, haplotype estimates obtained by PHASE v2 from all samples together will also
tend to be slightly more accurate, on average, than thoseobtained by analyzing the two population samples
separately.
Discussion
We have introduced a new algorithm for inferring haplotypes from population samples and have compared itin three ways with existing approaches. All comparisonspoint to the benefits of using models based on the coalescent to capture underlying patterns of haplotype
variation. Although one might expect the effects of taking the decay of LD into account to be greatest overlarger regions, the results shown in table 3 suggest that
there are gains to be made even at the level of datacollected over tens of kilobases. Further, the consistentimprovement in accuracy we observed across all comparisons argues that the use of these kinds of models
will be beneficial for most data sets—from most regionsof the genome—and not beneficial only “on average.”
Nevertheless, some users may worry that specific factors
for their data set (such as unusual patterns of polymorphism due to selection or complex marker-ascertainment
schemes) might make our approximate-coalescent model
inappropriate. We suspect this to be rare but note that
our experiments with missing data suggest a simple general strategy for comparing the likely accuracy of different haplotype-inference methods in such cases: deletea small proportion (e.g., 5%–10%) of the genotype data

Table 3Comparison of Accuracy of Methods for Imputing Missing Data
MISSINGNESS PATTERN FOR
Alleles Analyzed
Independently

Whole Genotypes

ALGORITHM

Separate

Combined

Separate

Combined

PHASE v2 recom (-MR)
PHASE v2 hybrid (-MQ)
PHASE v2 no recom (-MS)
Bayes-Dirichlet
PL-EM
HAP
hap2
Straw mansnphap

.033
.034
.049
.053
.086
NA
NA
.100
NA

.027
.031
.044
.049
.072
NA
NA
.102
NA

.044
.049
.066
.077
.088
.082
.101
.155
.188

.038
.041
.062
.063
.062
.085
.117
.158
.167

NOTE.—Each number is the genotype-imputation error rate (described intext). The “Separate” columns give the results for analysis of the 23 ED and
24 AD individuals separately; the “Combined” columns give the results foranalysis of all 47 individuals together. The differences between our new methods(PHASE v2-MR and -MQ) versus all other methods are statistically significant.
For the “Alleles Analyzed Independently” pattern, the results are based on a
total of 19,021 missing alleles; for the “Whole Genotypes” pattern, the results
are based on 18,958 missing alleles.

Stephens and Scheet: Accounting for Decay of Linkage Disequilibrium

and see which inference method produces more-accurate
genotype estimates (on average, over several repetitionsof the process). This strategy depends, of course, on theuse of a haplotype-inference method to impute missing
genotypes that is similar to the one they use to estimate
haplotypes. This is the case for all the methods implemented in PHASE—since they use the same statisticalmodel for both tasks—but may not be the case for allother methods. Thus, although, for PHASE methods, theresults in table 3 should reflect their relative accuracy inboth imputing missing genotypes and inferring haplotypes, other methods may reflect relative accuracy onlyin the specific task of imputing missing genotypes.
The high accuracy with which our method can impute
missing genotypes suggests that it could play a usefulrole in selecting which SNPs to genotype in large-scaleassociation studies. As pointed out by Johnson et al.
(2001), genotyping costs in such studies can be substantially reduced—without much loss in power—by
genotyping only an appropriately chosen subset of SNPs
(usually referred to as “tagSNPs”). One way of approaching this problem is to choose the tagSNPs so thatgenotypes at all other SNPs can be predicted with a highdegree of accuracy from the typed SNPs. Most published methods along these lines involve fairly simple
approaches to measuring how accurately genotypes at
a subset of SNPs predict the genotypes at others (e.g.,
Chapman et al. [2003] use linear regression). It would
be interesting to compare the accuracy of these simple
approaches for predicting missing genotypes with ourapproach. If our approach provides more-accurate predictions, then, when coupled with appropriate analyses
for identifying associations between genotype and phenotype, this could lead to more-efficient studies that
require fewer tagSNPs.
For the data considered here, genotype-imputation
estimates were more accurate for analysis of AA andED samples together than for separate analysis. This isconsistent with previous observations that haplotype inference methods are not especially sensitive to deviations from the assumption of Hardy-Weinberg equilibrium that underlies most of them (Fallin and Schork
2000; Stephens et al. 2001). However, in some studies,
haplotype estimates may tend to be more accurate if
samples from individuals of different ethnic backgrounds
are analyzed separately rather than together. Specifically, this might be expected to occur if the different
ethnicities are highly diverged in the region of study
and/or the sample sizes available for each ethnicity arelarge. As discussed above, randomly deleting a fractionof known genotypes and determining which analysis
approach (together or separate) provides more-accurate
genotype estimates provides a strategy for deciding
which approach is likely to provide the most-accurate
haplotype estimates in a particular case. Note, though,459

that this strategy addresses only the expected accuracyof haplotype estimates and not issues of bias: we would
always expect that analyzing samples from different
backgrounds together would tend to systematically underestimate differences in haplotype frequencies among
the groups, whereas analyzing samples separately would
tend to systematically overestimate these differences.
Although the algorithms implemented in PHASE v2
are among the most-accurate algorithms for haplotype
estimation available, they are also among the most computer intensive. The new algorithm we present here is
the most computer-intensive of all of these. The mostobvious reason for this is that computation of the conditional distribution p in expression (1) is substantiallymore time consuming than analogous computations
that ignore recombination. However, another importantfactor is that allowance for recombination tends to increase the number of “plausible” haplotypes for eachindividual, which also increases the computational cost.
Given the similar accuracy achieved by the other PHASE
methods, particularly on the simulated X chromosome
data, one might wonder whether the additional computing expenditure is worthwhile. Our view is that one
should use the method that provides the best results
within the time frame available. A single run of theslower of our new algorithms (-MR) on the largest of
the data sets that we considered here (the data fromDaly et al. [2001], with 103 SNPs for 129 individuals)
took roughly 3 h of computing time on a 3-GHz CPU
desktop machine, which remains well within the boundsof what we consider reasonable.
Results for the X chromosome data in table 2 highlight the fact that, in some settings, haplotypes estimatedby even the best available methods may be wrong forthe majority of sampled individuals. The low overall
accuracy in this case is presumably due to several factors, including the small sample size, the large numberof markers, and the fact that the markers are spreadover reasonably large regions. These results might lead
one to question the relevance of statistical methods for
estimating haplotypes in such settings. Fortunately,
however, for many analyses, conclusions will not be
critically dependent on every detail of each individual’s
haplotypes. (Indeed, it seems prudent to design analysesof statistically reconstructed haplotypes in such a waythat this is the case.) For this reason, although measuring the accuracy—by various criteria—of estimatedhaplotypes provides a convenient way to compare therelative merits of different algorithms, it seems impossible to translate such results into an absolute measureof their usefulness. In practice, what is needed is a way
to determine whether, for a particular data set, haplotypes can be estimated with sufficient accuracy to draw
reliable conclusions. For this, it seems essential that haplotype-reconstruction algorithms provide not only point

460

estimates for the haplotypes but also a reliable measureof the uncertainty in estimated haplotypes. As noted by
Stephens and Donnelly (2003), Bayesian approaches
seem to have an advantage here, since the Bayesianframework provides a natural way of assessing uncertainty, via the conditional distribution of the haplotypes,
given the genotype data.
Among Bayesian approaches, the new algorithms we
present here should better capture the uncertainty inestimated haplotypes than do previous methods, because
the assumptions underlying our approach are more realistic (as evidenced by the reduced error rates shownin tables 2 and 3). In particular, explicit allowance forrecombination tends to increase the number of “plausible” haplotypes for an individual and thus leads togreater uncertainty in the estimated haplotypes. It may
seem counterintuitive to claim that a method withgreater uncertainty in estimated haplotypes is superior
to a method with less uncertainty. However, methods
that provide for less uncertainty risk being overconfident in their assessments of estimated haplotypes, withthe resulting danger of being overconfident in subsequent conclusions. There are, however, situations forwhich our new algorithm, as we have applied it here,

Am. J. Hum. Genet. 76:449–462, 2005

may be underconfident in its estimated haplotypes. Thiscan occur when analyzing small data sets (few markers
and/or few individuals), because, for such data sets, themethod may tend to overestimate the recombination
rates: with few markers, inference of the recombination
rate will depend heavily on the prior distribution for
the recombination parameter, and the rather flat priorthat we use here perhaps allows too much weight on
unrealistically large values. Fortunately, this problem iseasily solved by use of a more realistic prior for therecombination rate, although appropriate choice of thisprior will typically depend on the organism and theregion under study; we wanted to avoid this additional
subjectivity here. The prior can, however, be altered asa parameter in the PHASE v2 software package, whichis available from the Stephens Lab Web site.
Acknowledgments
We thank two anonymous reviewers for their helpful andconstructive comments. This work was supported by National Institutes of Health grants 1RO1HG/LM02585-01(to M.S.) and T32 HG00035 (to P.S.).

Appendix A
Further Algorithmic Details
We refer readers to the publication by Stephens and Donnelly (2003) for details of the PL and for more detailsof that aspect of the algorithm. We focus on providing further details of steps 1 and 2 of our algorithm, outlinedin the text.
Step 1: Updating Individuals
Updating the haplotype pair of individual i (including alleles at loci with missing genotypes) proceeds as follows:
1. Make a list Hi of all haplotype pairs that are both compatible with the observed genotypes of individual i,and occur in the current list L of plausible haplotypes (this list is created by the PL procedure; see Stephens andDonnelly [2003] for more details).
2. For each haplotype pair in Hi, compute the probability of the pair using expression (1).
3. Sample a random pair in Hi according to these probabilities and set the haplotypes of i to this pair.
4. Impute missing positions: for each of the two haplotypes of i (one at a time), sample alleles at the positions
with missing genotypes from their conditional distribution, given the current imputed values of the alleles at positions
with nonmissing genotypes. This can be done efficiently, by using the forward-backward algorithm (Rabiner 1989)
to first simulate the values of the hidden Markov chain (X,T) for the haplotype and then imputing the missing
alleles independently at each locus, conditional on the simulated value of (X,T) . If the new haplotype pair of i is
not in Hi, add it.
Note that step 4 is not strictly necessary, in that the missing alleles are also imputed by the other steps. It was
included to improve mixing (although we have not performed a detailed study of its effectiveness).

461

Stephens and Scheet: Accounting for Decay of Linkage Disequilibrium

Step 2: Updating Recombination Parameters
The recombination parameters, r̄ and l p (l 1 , … ,lL⫺1) , are initialized to be r̄ p 0.0004 (on the basis of approximate genomewide average for humans) and l l p 1. They are then updated, using Langevin MH proposals
on a log scale (Besag 1994). Langevin MH updates are MH updates that use information from the derivative of
the posterior density to propose moves in “sensible” directions that are likely to have a reasonable chance of being
accepted. The standard MH acceptance probability is used. Specifically, the update proposals for r̄ and l are:
1. Update r̄ by proposing to add 0.5jr¯2d ⫹ e to log (r̄), where e ∼ N(0,jr¯2) and d p 冘l (r̄l l D l) and D l is the partial
derivative of the PAC-B likelihood (Li and Stephens 2003), with respect to rl p l lr̄ . (These derivatives can be
efficiently computed analytically by use of the forward–backward algorithm for hidden Markov models (Rabiner1989).
2. Update l, by proposing to add 0.5jl2dl ⫹ e to each log (l l), wheredl p r̄l l D l ⫺log (l l)
1.15 2(the last term coming from the N(0,1.15 2 ) prior on log (l l)) and e ∼ N(0,jl2 ).
Here, jr¯ and jl control the expected size of proposed jumps and need to be tuned to the particular data set being
analyzed, to ensure reasonable mixing behavior. We adopt a rather simplistic procedure to perform this automatictuning that is based on the idea that jr¯ and jl should be set so that acceptance probabilities are not too close to
0 or 1. Initially, we set jr¯ p jl p 1. Then, halfway through the burn-in iterations, better values can be found by
repeating the following until acceptance rates for both r and l are in the range 0.3–0.7:
1. update r̄,l 10 times each, using current values of jr¯ and jl, and
2. if the proportion of acceptances of r̄ (respectively l) was not in the range 0.3–0.7, then divide or multiplyjr¯ (respectively jl) by 1 ⫹ u, where u is random uniform on [0,1).
If necessary, this tuning can be repeated after the burn-in iterations have been completed, before the main iterations
are performed.

Appendix B
List of Genes Used for Results Shown in Table 3
Genes are listed by HUGO Gene Nomenclature Committee abbreviation: CD36, CEBPB, CRF, CRP, CSF3,
CYP4A11, CYP4F2, DCN, EPHB6, F11, F2RL3, F3, IFNG, IGF2, IL1B, IL2RB, IL4, IL5, IL10A, IL11, IL19,
IL20, IL21R, IL22, IL24, KEL, LTA, LTB, MAP3K8, MMP3, PFC, PLAUR, PLG, PROC, PTGS2, SCYA2, SELE,
SELL, SELPLG, SERPINC1, SFTPA1, SMP1, STAT6, TFPI, THBD, TNF, TNFAIP2, TNFAIP3, TRAF6, andVCAM1.

Electronic-Database Information
Accession numbers and URLs for data presented herein are
as follows:

David Clayton Software page, http://www-gene.cimr.cam.ac.uk/clayton/software/
HAP Webserver,http://www1.cs.columbia.edu/compbio/hap
HUGO Gene Nomenclature Committee, http://www.gene.ucl.ac.uk/nomenclature/
Stephens Lab Web site, http://www.stat.washington.edu/
stephens/software.html (for PHASE)

UW-FHCRC Variation Discovery Resource, http://pga.gs.washington.edu

References
Besag JE (1994) Discussion on the paper by Grenander andMiller. J R Stat Soc B 56:591–592
Carlson C, Eberle M, Rieder M, Smith J, Kruglyak L, Nickerson D (2003) Additional SNPs and linkage-disequilibrium
analyses are necessary for whole-genome association studies
in humans. Nat Genet 33:518–521
Chapman J, Cooper J, Todd J, Clayton D (2003) Detecting
disease associations due to linkage disequilibrium using hap-462
lotype tags: a class of tests and the determinants of statisticalpower. Hum Hered 56:18–31
Clark AG (1990) Inference of haplotypes from PCR-amplified
samples of diploid populations. Mol Biol Evol 7:111–122
Crawford D, Bhangale T, Li N, Rieder M, Nickerson D, Stephens M (2004) Evidence for substantial fine-scale variationin recombination rates across the human genome. Nat Genet
36:700–706
Daly M, Rioux J, Schaffner S, Hudson T, Lander E (2001) Highresolution haplotype structure in the human genome. Nat
Genet 29:229–232
Eskin E, Halperin E, Karp R (2003a) Efficient reconstructionof haplotype structure via perfect phylogeny. J Bioinform
Comput Biol 1:1–20——— (2003b). Large scale reconstruction of haplotypes fromgenotype data. In: Proceedings of the seventh annual international conference on research in computational molecularbiology (RECOMB 2003), pp 104–113
Evans G (1993) Practical numerical integration. Wiley and Sons,
Chichester, United Kingdom
Excoffier L, Laval G, Balding DN (2003) Gametic phase estimation over large genomic regions using an adaptive window
approach. Hum Genomics 1:7–19
Excoffier L, Slatkin M (1995) Maximum-likelihood estimationof molecular haplotype frequencies in a diploid population.
Mol Biol Evol 12:921–927
Fallin D, Schork NJ (2000) Accuracy of haplotype frequency
estimation for biallelic loci, via the expectation-maximization algorithm for unphased diploid genotype data. Am J
Hum Genet 67:947–959
Fearnhead PN, Donnelly P (2001) Estimating recombination
rates from population genetic data. Genetics 159:1299–1318
Greenspan G, Geiger D (2004) Model-based inference of haplotype block variation. J Comput Biol 11:493–504
Halperin E, Eskin E (2004) Haplotype reconstruction fromgenotype data using imperfect phylogeny. Bioinformatics 20:
1842–1849
Hawley M, Kidd K (1995) HAPLO: a program using the EM
algorithm to estimate the frequencies of multi-site haplotypes.
J Hered 86:409–411
Heckerman D, Chickering DM, Meek C, Rounthwaite R, Kadie C (2000) Dependency networks for inference, collaborative filtering, and data visualization. J Machine Learning
Res 1:49–75
Hudson RR (1991) Gene genealogies and the coalescent process. In Futuyma D, Antonovics J (eds) Oxford surveys inevolutionary biology. Vol 7, pp 1–44. Oxford University
Press, Oxford, United Kingdom

Am. J. Hum. Genet. 76:449–462, 2005

Johnson GC, Esposito L, Barratt BJ, Smith AN, Heward J,
Genova GD, Ueda H, Cordell HJ, Eaves IA, Dudbridge F,
Twells RC, Payne F, Hughes W, Nutland S, Stevens H, Carr
P, Tuomilehto-Wolf E, Tuomilehto J, Gough SC, Clayton DG,
Todd JA (2001) Haplotype tagging for the identification ofcommon disease genes. Nat Genet 29:233–237
Kimmel G, Shamir R (2004) Maximum likelihood resolutionof multi-block genotypes. In: Proceedings of the Eighth Annual International Conference on Research in Computational Molecular Biology (RECOMB’04). The Association
for Computing Machinery, San Diego, pp 2–9
Li N, Stephens M (2003) Modelling linkage disequilibrium,
and identifying recombination hotspots using SNP data. Genetics 165:2213–2233
Lin S, Chakravarti A, Cutler DJ (2004) Haplotype and missing
data inference in nuclear families. Genome Res 14:1624–1632
Lin S, Cutler DJ, Zwick ME, Chakravarti A (2002) Haplotype
inference in random population samples. Am J Hum Genet
71:1129–1137
Long JC, Williams RC, Urbanek M (1995) An E-M algorithm
and testing strategy for multiple locus haplotypes. Am J Hum
Genet 56:799–810
McVean GAT, Myers SR, Hunt S, Deloukas P, Bentley DR,
Donnelly P (2004) The fine-scale structure of recombination
rate variation in the human genome. Science 304:581–584
Niu T, Qin ZS, Xu X, Liu JS (2002). Bayesian haplotype inference for multiple linked single-nucleotide polymorphisms.
Am J Hum Genet 70:157–169
Ptak SE, Roeder AD, Stephens M, Gilad Y, Pääbo S, Przeworski M (2004) Absence of the TAP2 human recombination hotspot in chimpanzees. PLoS Biol 2:849–855
Qin ZS, Niu T, Liu JS (2002) Partial-ligation–expectationmaximization algorithm for haplotype inference with singlenucleotide polymorphisms. Am J Hum Genet 71:1242–1247
Rabiner LR (1989) A tutorial on hidden Markov models andselected applications in speech recognition. Proc IEEE 77:
257–286
Schaffner SF (2004) The X chromosome in population genetics. Nat Rev Genet 5:43–51
Stephens M, Donnelly P (2000) Inference in molecular population genetics. J R Stat Soc Ser B 62:605–655
Stephens M, Donnelly P (2003) A comparison of Bayesian
methods for haplotype reconstruction from population genotype data. Am J Hum Genet 73:1162–1169
Stephens M, Smith NJ, Donnelly P (2001) A new statisticalmethod for haplotype reconstruction from population data.
Am J Hum Genet 68:978–989