Identifying nineteenth century genealogical links from genotypes

Abstract

We have developed a likelihood method to identify moderately distant genealogical relationships from genomewide scan data. The aim is to compare the genotypes of many pairs of people and identify those pairs most likely to be related to one another. We have tested the algorithm using the genotypes of 170 Tasmanians with multiple sclerosis recruited into a haplotype association study. It is estimated from genealogical records that approximately 65% of Tasmania’s current population of 470,000 are direct descendants of the 13,000 female founders living in this island state of Australia in the mid-nineteenth century. All cases and four to five relatives of each case have been genotyped with microsatellite markers at a genomewide average density of 4 cM. Previous genealogical research has identified 51 pairwise relationships linking 56 of the 170 cases. Testing the likelihood calculation on these known relative pairs, we have good power to identify relationships up to degree eight (e.g. third cousins once removed). Applying the algorithm to all other pairs of cases, we have identified a further 61 putative relative pairs, with an estimated false discovery rate of 10%. The power to identify genealogical links should increase when the new, denser sets of SNP markers are used. Except in populations where there is a searchable electronic database containing virtually all genealogical links in the past six generations, the algorithm should be a useful aid for genealogists working on gene-mapping projects, both linkage studies and association studies.

Appendix: Further details of the approximate likelihood calculation

Consider a pair of degree d cousins and the corresponding pair of related haplotypes (inherited from the cousins’ related parents) for a particular chromosome. Denote the two haplotypes h₁=(h₁₁,h₁₂,...) and h₂=(h₂₁,h₂₂,...), where h _ik is the allele at marker k along the ith haplotype . Then to evaluate the numerator in Eq. 2 we need to evaluate $\Pr((h_1,h_2)|d).$ We consider integer values of d between 3 (first cousins) and 13 (e.g. 6th cousins) inclusive.

Define the IBD process {D^d} for h₁ and h₂ by

$$ D^d_k = \left\{ \begin{array}{*{20}l} 1 & \hbox{if } h_1 \hbox{ and } h_2 \hbox{ are IBD at marker } k \\ 0 & \hbox{if } h_1 \hbox{ and } h_2 \hbox{ are not IBD at marker } k \end{array} \right.$$

Denote the population frequencies of the alleles h_1k and h_2k by f_1k and f_2k respectively. Suppose that there is a probability ɛ that one of the alleles h_1k and h_2k has been incorrectly genotyped or inferred, or that one of the alleles has undergone a mutation in the meioses separating the pair of relatives. Then, we have the following probabilities for the genotypes at marker k, given the IBD status at marker k:

$$\Pr((h_{1k},h_{2k})|D^d_k=0) = f_{1k} f_{2k}$$

(5)

$$\Pr((h_{1k},h_{2k})|D^d_k=1) = \left\{ \begin{array}{*{20}l}\varepsilon f_{1k} f_{2k,} & h_{1k} \neq h_{2k} \\ (1-\varepsilon) f_{1k} + \varepsilon f_{1k}^2, & h_{1k} = h_{2k}.\end{array}\right. $$

(6)

When the allele on a particular haplotype at a particular marker is ambiguous, the probabilities Eqs. 5 and 6 were evaluated using both equally likely alleles and averaged.

The IBD process {D^d} is not Markov for the cousin relationships we consider. However, it is possible to construct an augmented process {A^d} which contains all of the information of the IBD process {D^d} and which is Markov in the absence of interference (Donnelly 1983). McPeek and Sun (2000) present the minimal augmented Markov process {A³} for first cousins in detail: they define 7 states 1,...,7, with states A³=5 and A³=7 corresponding to D³=1. More generally, suppose that the haplotypes h₁ and h₂ carried by dth degree cousins are descended from a pair of related first cousin haplotypes g₁ and g₂ via an additional d−3 meioses. Then the probability of the haplotype h₁ being IBD with the haplotype g₁ and the haplotype h₂ being IBD with the haplotype g₂ at a particular marker is 1/2^d-3: all d−3 of these additional meioses must be “switched” a certain way. Define the Markov process {S^d-3}, which counts how many of these meioses are switched the right way for haplotypes h₁ and h₂ to be IBD with haplotypes g₁ and g₂, respectively. {S^d-3} has d−2 possible states: 0,1,...,d−3. Then the product of {S^d-3} and the process {A³} for the first cousin haplotypes g₁ and g₂ is an augmented Markov process {A^d} which contains all of the information of the non-Markov IBD process {D^d}.

The augmented process {A^d} has 7(d−2) states 1,...,7(d−2). To move from marker k−1 to marker k in the likelihood calculation requires evaluation of [7(d−2)]² transition probabilities

$$ \Pr(A^{d}_k=j | A^{d}_{k-1}=i, \theta_{k-1,k}), \quad i,j = 1,\ldots,7(d-2) $$

where θ_k-1,k is the recombination fraction between markers k−1 and k. To reduce computation we follow the suggestion of McPeek and Sun and approximate the IBD process {D^d} by a Markov process {B^d} with the correct transition probabilities for adjacent markers: i.e. set

$$ \Pr(B^{d}_k=j | B^{d}_{k-1}=i) = \Pr(D^{d}_k=j | D^{d}_{k-1}=i), \quad i,j = 0,1. $$

First, evaluate

$$\Pr(D_k^d=1 | D_{k-1}^d=1) = \Pr(D_k^3 =1 | D_{k-1}^3=1) \Pr(S_k^{d-3}=d-3 | S_{k-1}^{d-3}=d-3).$$

(7)

Now, for the process {S^d-3} to stay in state d−3 with all d−3 meioses switched a certain way, there must be no recombinations between markers k−1 and k. Thus

$$\Pr(S_k^{d-3}=d-3 | S_{k-1}^{d-3}=d-3) = (1-\theta)^{d-3},$$

(8)

where θ=θ_k-1,k.To evaluate the probability that the process {D³} stays in state 1 (with the first cousin haplotypes staying IBD) between markers k−1 and k, refer back to the augmented Markov process {A³}, which contains all the information of {D³}. Using the labelling of the seven states of {A³} and the transition probabilities between them given by McPeek and Sun (2000), we have

$$ \begin{aligned} \Pr(D_k^3 =1 | D_{k-1}^3=1) &= \sum_{i,j \in \{A^3: D^3=1\}} \Pr(A_{k-1}^3=i | D_{k-1}^3=1) \Pr(A_k^3=j | A_{k-1}^3=i) \\ &= \frac{{1}}{{2}}\left[\Pr(A_k^3=5 | \Pr(A_{k-1}^3=5) + \Pr(A_k^3=7 | \Pr(A_{k-1}^3=5) \right. \\ & \quad + \left. \Pr(A_k^3=5 | \Pr(A_{k-1}^3=7) + \Pr(A_k^3=7 | \Pr(A_{k-1}^3=7) \right] \\ &= \frac{{1}} {{2}}\left[(1-\theta)^2 \psi^2 + \theta^2 \phi^2 + 2 \psi^2 \phi + \psi^3 \right],\end{aligned}$$

(9)

where ψ=θ²+(1−θ)² and ϕ=1 − ψ. Substituting Eqs. 8 and 9 into Eq. 7 gives

$$\Pr(D_k^d=1 | D_{k-1}^d=1) = \frac{{1}}{{2}} (1-\theta)^{d-3} \left[(1-\theta)^2 \psi^2 + \theta^2 \phi^2 + 2 \psi^2 \phi + \psi^3 \right].$$

(10)

For convenience, write the transition probability from state i to j in abbreviated notation: $t_{ij} = \Pr(D_k^d=j | D_{k-1}^d=i)$. Equation 10 gives t₁₁. Then

$$ t_{10}=1-t_{11}.$$

Now the Markov chains {A³} and {S^d-3} have stationary distributions

$$\pi_1 = \pi_2 =\pi_4 =\pi_5 =\pi_6 =\pi_7 = 1/8, \ \pi_3=1/4$$

and

$$\pi_{\ell} =\frac{{(d-3)!}}{{\ell! (d-3-\ell)!}}\frac{{{1}}}{{2^{d-3}}}, \quad \ell=0,1,\ldots,d-3$$

respectively. Thus the process {D^d}, which has value 1 when S^d-3=d−3 and A³=5 or A³=7, has stationary distribution

$$\pi_0 = 1-\frac{{1}}{{2^{d-1}}} {,\ } \pi_1 = \frac{{1}}{{2^{d-1}}} $$

(11)

satisfying

$$\left( \begin{array}{*{20}l} t_{11} & t_{01} \cr t_{10} & t_{00} \cr \end{array}\right) \left( \begin{array}{*{20}l} \pi_1 \cr \pi_0 \end{array}\right)=\left( \begin{array}{*{20}l} \pi_1 \cr \pi_0 \end{array} \right).$$

(12)

Substituting Eq. 11 into Eq. 12 and solving for t₀₁ gives

$$ t_{01} = {{1 - t_{11}}\over {2^{d-1}-1}}, $$

and, finally,

$$ t_{00} = 1 - t_{01}. $$