Identification of candidate IBD segments

When applying the GERMLINE algorithm to detect candidate IBD segments, we do not permit any mismatching alleles in the shared haplotype. Each candidate IBD segment is defined by its starting and ending genome coordinates, the pair of sample identifiers, and the haplotype index (1 or 2) of the shared haplotype for each sample.

The ibdwindow parameter in Beagle version 4 determines the number of markers included in each window when using the GERMLINE algorithm to find candidate IBD segments. The ibdwindow parameter is equivalent to the GERMLINE bits parameter. Too large a value may result in missing short segments of IBD, while too small a value will increase computation time. The default value of 64 is suitable for SNP arrays with 1 million SNPs across the genome, as at this marker density, 64 markers correspond to ~0.2 cM, which is significantly shorter than the default threshold on IBD segment length. For SNP array data, we recommend setting this parameter to approximately the average number of markers per 0.2 cM.

The ibdcm parameter in Beagle version 4 controls the minimum genetic length of a candidate IBD segment. A value that is too small will result in increased computing time while not contributing much to IBD detection as small candidate segments are unlikely to pass the LOD score threshold. The default value of 1.0 cM was chosen based on the relatively low power to detect smaller segments in SNP array data (see Results).

Haplotype frequency models for IBD and non-IBD

We start with a model for haplotype frequencies. We use the Beagle HMM (S. R. Browning and B. L. Browning 2007), but our approach is general and it could be readily adapted to other HMMs for haplotype frequencies. The HMM for haplotype frequencies determines a HMM for unrelated individuals and a HMM for parent–offspring pairs (Browning and Browning 2009). We calculate the probability of the observed genotype data for a pair of individuals under a non-IBD model (the likelihood of non-IBD), using the HMM for unrelated individuals, and we calculate the probability of the observed genotype data for a pair of individuals under an IBD model (the likelihood of IBD), using the HMM for a parent–offspring pair since a parent–offspring pair also shares one haplotype identical by descent.

A HMM is defined by its state space, initial probabilities, transition probabilities, set of emitted symbols, and emission probabilities (Rabiner 1989). Since the Beagle HMMs for unrelated individuals and parent–offspring pairs have been fully described previously (S. R. Browning and B. L. Browning 2007, 2009), we give only a brief description of these HMMs here. In the Beagle HMM, there is a set of haploid hidden states S_m corresponding to each marker m. Each haploid hidden state corresponds to a cluster of haplotypes that are locally similar around marker m. The hidden state for a diploid individual at marker m is an element $(s_1,s_2)\in S_m^2$, where s₁ is the state of the first haplotype, and s₂ is the state of the second haplotype of the individual. The hidden state of a parent–offspring pair at marker m is an element $(s_1,s_2,s_3)\in S_m^3$, where (s₁, s₂) is the hidden state of the parent, (s₁, s₃) is the hidden state of the offspring, and s₁ is the hidden state of the shared haplotype.

In the Beagle HMM, each haploid hidden state at marker m is labeled with one of the marker’s alleles (more than one hidden state at marker m can be labeled with the same allele). The emission probability for the labeled allele is 1, while the emission probability for any other allele is 0. For parent–offspring pairs, this means that the probability of the observed genotypes given the hidden state $(s_1,s_2,s_3)\in S_m^3$ is 1 if the alleles labeling s₁ and s₂ are consistent with the observed genotype of the parent at marker m and the alleles labeling s₁ and s₃ are consistent with the observed genotype of the offspring at marker m; the probability is 0 otherwise.

Each haplotype used to build the Beagle HMM has a unique path through the model. Thus each haploid state has an associated count of how many haplotypes pass through the state, and similarly each possible transition has an associated count. Consider a haploid transition in the Beagle HMM. Call the state that the transition starts from at marker m the “source” state and the state that the transition goes to at marker m + 1 the “destination” state. The transition probability is the count associated with that transition divided by the count associated with the source state for the transition. In other words, of those haplotypes passing through the source state, the transition probability is the proportion that transitions into the destination state. Given a starting marker m, the initial probability for a haploid state at marker m is equal to the proportion of haplotypes that pass through that state, which is the count associated with that state divided by the total number of haplotypes used to build the model. Transition and initial probabilities for unrelated diploid individuals (pairs of hidden states) or parent–offspring pairs (triples of hidden states) are obtained by multiplying the corresponding haploid probabilities.

Since IBD segments are typically much shorter than their corresponding chromosome, we reduce computation time by calculating likelihoods using only genotype data in the interior of the candidate IBD segment. In this way, we can also avoid modeling the recent recombination events that demarcate the boundaries of the IBD segment. It is difficult to identify the IBD segment endpoints with high accuracy. Incorrectly including some non-IBD markers at the ends of a real IBD segment can result in a severely reduced LOD score for the segment, while incorrectly removing small parts of the ends of a real IBD segment tends to result in only a small drop in LOD score. We thus trim a small fixed number of markers from each end of the candidate IBD segment and compute likelihoods of the non-IBD and IBD models in the trimmed genomic interval. The trimmed markers are restored to the IBD segment after calculation of the likelihoods. The number of markers to trim when calculating the LOD score is controlled by the ibdtrim parameter in Beagle version 4. Increasing the trim number will reduce power to detect short segments. When a short candidate segment is trimmed there may not be enough markers left to provide sufficient information to be confident of IBD. The default trim value (40) was chosen based on analyses of the Wellcome Trust Case Control Consortium 2 data, by looking for a value that would maximize the amount of IBD detected (data not shown). For SNP array data we recommend setting the trim value to approximately the average number of markers per 0.15 cM.

The likelihoods for the IBD and non-IBD models are calculated using Baum’s forward algorithm (Baum 1972). For the non-IBD model, the probability of the observed data for a pair of individuals is the product of the probabilities for each individual.

Other improvements to Beagle in version 4

As well as implementing Refined IBD, Beagle version 4 has improvements to haplotype phasing and to general usability. Beagle version 4 reports the consensus haplotype for each individual. Other haplotype phasing programs also use consensus haplotypes (Scheet and Stephens 2006; Li et al. 2010). Previous versions of Beagle have used the Viterbi algorithm (Viterbi 1967) to generate the reported phased haplotypes. However, we have found the consensus haplotypes to be more accurate than the haplotypes obtained from the Viterbi algorithm.

For consensus haplotypes, haplotypes are estimated at multiple iterations (or from multiple runs) of the phasing algorithm, and the results are merged. The top row of Figure 1 illustrates the procedure for obtaining consensus haplotypes in Beagle version 4. The haplotype phasing module involves multiple iterations of estimating (sampling) haplotypes based on a provisional model and then updating the model based on the new estimated haplotypes. By default, four pairs of haplotypes are sampled per individual per iteration. Haplotypes estimated in the first few iterations are not likely to be very accurate, because the provisional model is still in the initial stages of converging toward a good solution. Thus, in Beagle version 4, the consensus haplotypes are obtained from all sampled haplotypes after a specified number of burn-in iterations (five burn-in iterations and five additional iterations by default). Under default settings there are 20 pairs of sampled haplotypes per individual that are used for obtaining the consensus haplotypes (4 pairs per iteration × 5 iterations after burn-in).

The first step in obtaining the consensus haplotypes is to obtain consensus genotypes for those genotypes that were missing. Consensus genotypes are obtained by taking the most frequently sampled genotype, breaking ties randomly. After consensus genotypes are obtained, the consensus phasing for an individual is obtained by working along the chromosome, one pair of successive heterozygous genotypes at a time. A pair of successive heterozygous genotypes has no intervening heterozygous genotypes. Only sampled haplotype pairs having the consensus heterozygous genotype at both markers are used to determine the consensus phasing. The consensus phasing of two successive heterozygous genotypes is determined by majority vote, breaking ties randomly. For example, in phasing the genotypes AC and TG, if 16 sampled haplotype pairs have AT/CG phase, while 4 sampled haplotype pairs have AG/CT phase, the consensus haplotypes will have the AT/CG phase.

Beagle version 4 uses Variant Call Format for input and output data files (Danecek et al. 2011). Variant Call Format is a standard, widely used format for genotype data (1000 Genomes Consortium 2010), and use of this format will reduce the need for tedious data file format conversion.

Beagle version 4 also uses a sliding marker window that makes the memory usage independent of the number of markers in the data set. Decreased accuracy near the edge of the marker windows is avoided by using overlapping windows and trimming half of the overlap from each window prior to merging data in adjacent windows. Haplotypes in adjacent marker windows are aligned using a heterozygote near the middle of the overlap.

Scale factors

The Beagle HMM is represented by a directed acyclic graph. When the model is constructed, a process of node merging occurs. Two nodes of the graph, x and y, are merged if the maximum difference in downstream frequencies is less than a threshold (B. L. Browning and S. R. Browning 2007). The threshold is

where m is the scale factor, b is the shift parameter, and n_x and n_y are the numbers of haplotypes whose path through the graph includes nodes x and y, respectively. The scale and shift parameters were originally introduced to control the degree of parsimony of the fitted Beagle model when performing association testing (B. L. Browning and S. R. Browning 2007). Larger values of these parameters result in more merging of nodes and hence a more parsimonious model. Without the added parsimony, each haplotype cluster may contain very few observations, reducing power to detect an association.

For Beagle’s other applications (phasing, imputation, and IBD detection), we have not found the shift parameter to be useful, so we assume a shift parameter of b = 0 for the remainder of this section. For phasing and imputation, we have found that a scale factor of m = 1 performs well. For IBD detection, a different scale factor (the IBD scale factor) can be used for the final model while continuing to use a scale factor of 1 for the haplotype phasing step. We have previously used IBD scale factors of 1 (Browning and Browning 2010) and 2 (B. L. Browning and S. R. Browning 2011). Recently we have realized that for a given sample size (number of genotyped individuals), the choice of IBD scale parameter for values ≥2 is somewhat arbitrary, provided that appropriate compensatory adjustment is made to the threshold for significance of the fastIBD score or the Refined IBD LOD score (data not shown). However, as the sample size changes, the optimal choice of IBD scale factor for a fixed score threshold changes. For a given choice of IBD scale parameter, as more individuals are added to the data, the fitted model becomes larger (less parsimonious). A larger model allows for higher precision in haplotype frequency estimation, resulting in fewer false positives. On the other hand, the requirement that shared haplotypes must traverse the same path through the model to be declared IBD becomes more onerous as the model size increases, and detection rates can drop.

Therefore, to have a single LOD score threshold regardless of sample size, the IBD scale factor must increase as the sample size increases and the size (complexity) of the fitted Beagle model must stay approximately constant, to maintain power to detect IBD as sample size increases. Since model complexity is controlled by the merging threshold given in Equation 1, when the sample size increases by a factor of k, the IBD scale factor can be increased by a factor of $\sqrt{k}$ to keep the typical threshold at approximately the same level, resulting in a similar size of model. We chose to make the default setting for the IBD scale factor $\sqrt{n/100}$, when the sample size, n, is >400. This results in IBD scale factors of 2.2 for 500 samples, 4.5 for 2000 samples, and 7.1 for 5000 samples. For sample sizes <400, we set the default IBD scale factor to 2, as decreasing the IBD scale factor below 2 can reduce power to find IBD.

In summary, with reference to Equation 1 and Figure 1, a scale factor of 1 is used for model building during haplotype phasing (top row of Figure 1). However, an IBD scale factor >1 is used for building the final Beagle model for IBD refinement (bottom row of Figure 1), as described above. A shift of 0 is used in all cases.

Simulated SNP data

To assess false positive and true positive IBD detection rates, we simulated data. In our simulation data, we attempted to match both current and historical effective population sizes, to obtain a good match with real data. Ancient historical population sizes affect the number of common variants and the extent of LD between them. The extent of LD affects power to detect IBD and can affect false positive IBD detection rates. Current and very recent population sizes affect the number of rare variants and the amount of detectable IBD. The amount of detectable IBD in the simulation is critical. Our method is designed for large outbred populations in which the amount of detectable IBD is low, so we want the simulated data to reflect this.

Recent analyses of the allele frequency spectrum, with particular attention paid to the rare end of the spectrum, have shown that explosive population growth has occurred in the past few hundred generations (Keinan and Clark 2012). In our view, previous models fitted to sequence data do not go far enough in modeling this growth, as they allow only for a single rate of recent growth, resulting in growth rate estimates of 2% (Nelson et al. 2012) or 9% (Coventry et al. 2010) per generation. In contrast, census data show that the rate of population growth has accelerated in the past few hundred generations and is currently ~30% per generation globally (Keinan and Clark 2012).

We used Fastsimcoal (Excoffier and Foll 2011) to simulate sequence data that we then thinned to obtain simulated SNP array data. We simulated 10 regions, each with 30 Mb of sequence on 2000 diploid individuals. We used a mutation rate of 2.5 × 10⁻⁸ (Nachman and Crowell 2000) and a recombination rate of 10⁻⁸ (i.e., 1 Mb = 1 cM). Our simulation scheme was designed with European populations in mind, as our available real data are from European populations. The effective population size was initially (prior to expansion beginning 300 generations ago) 3000 diploid individuals. This reflects the European effective population size estimated using LD between common variants (Tenesa et al. 2007). In our simulations, the effective population size began to grow 300 generations ago (timing reflects the advent of large-scale organized agriculture) at a rate of 1.8% per generation (reflecting, e.g., the 1.7% growth rate estimate in Nelson et al. 2012), reaching 270,000 by 50 generations ago. We modeled population growth rate increases in the past 50 generations based on English census data, as shown in Supporting Information, Figure S1. At 50 generations ago, we increased the growth rate to 5%, giving effective population size 2 million at 10 generations ago. We increased the growth rate further, to 25% per generation, for the final 10 generations, yielding an effective population size of 24 million (2 × 10⁶ × exp(0.25 × 10)) at the current generation.

After generating the sequence data, we created simulated SNP array data from it by removing all variants with more than two alleles and all variants with frequency <2%, and by selecting variants from those remaining to obtain ~1000 variants per 30-Mb region (corresponding to a SNP density of 1 million SNPs genome-wide) with minor allele frequencies uniformly distributed between 2% and 50%. We then added genotype error at a rate of 0.05%, reflecting the very high accuracy seen in current genotyping arrays after applying standard quality control filters (Steemers et al. 2006). Genotype error was introduced by converting homozygote genotypes to heterozygote and by converting heterozygote genotypes to a randomly chosen homozygote genotype. We also removed haplotype phase information.

We performed analyses on all 2000 simulated individuals and on a subset of 500 individuals. We used Refined IBD with minimum segment length 0.5 cM and LOD score thresholds of 3 and 4, with the remaining parameters at their default settings. We ran fastIBD (Beagle version 3.3.1) (B. L. Browning and S. R. Browning 2011) with IBD scale at the default value of 2 for n = 500 (this is very close to our new recommendation of $\sqrt{n/100}=2.2$) and with IBD scale equal to $\sqrt{n/100}=4.5$ for n = 2000. We used fastIBD score thresholds of 10⁻⁸ and 10⁻¹⁰. We ran fastIBD with 10 different random-number seeds and merged the results. We ran Beagle IBD (Beagle version 3.3.1) (Browning and Browning 2010) with 10 different random-number seeds, with the default IBD scale of 2 for 500 individuals. We did not run Beagle IBD with all 2000 individuals due to its long computing times.

We ran GERMLINE version 1.5.1 (Gusev et al. 2009) with parameters used in Gusev et al. (2011). Specifically, we used options “-haploid -min_m 1 -bits 32 -err_hom 1 -err_het 1”. The “-min_m 1” option means that GERMLINE reports only IBD segments with estimated length ≥1 cM. The “-haploid” option ensures that GERMLINE makes use of the haplotype phase information. Our previous published analyses with GERMLINE used GERMLINE’s default setting that does not use haplotype phase information (B. L. Browning and S. R. Browning 2010, 2011). As seen by comparing the results presented here with the results in our earlier work, utilization of haplotype phase information greatly improves GERMLINE’s performance on SNP array data. We used the phased haplotypes output from the Beagle Refined IBD analysis as input to GERMLINE. These haplotypes are based on a consensus of haplotypes sampled at different iterations of the Beagle phasing algorithm, as described above, and have accuracy higher than that from previous versions of Beagle run with default settings. The high SNP density (equivalent to a SNP array with 1 million SNPs) also contributes to high phasing accuracy. The high accuracy of these haplotypes facilitates the strong performance of GERMLINE in these data.

We used the full simulated phase-known sequence data to determine the true IBD status, so that we could assess the accuracy of the IBD estimated from the thinned, phase-agnostic SNP data. When determining the true IBD, we ignored variants with ≤10 copies in the 2000 individuals, as very recent mutations disrupt sequence identity. For two IBD haplotypes separated by m meioses, looking along the chromosome, recombination ends the IBD at rate m times the recombination rate, while new mutations disrupt the sequence identity at rate m times the mutation rate. In our simulation, the mutation rate is 2.5 times the recombination rate, so in any IBD segment we expect an average of 2.5 identity-disrupting new mutations. Ignoring variants with ≤10 copies, we declared pairs of haplotypes with identical sequence for at least 0.1 cM to be segments of true IBD for the purpose of assessing IBD detection accuracy.

Simulated sequence data

Using the simulated phased sequence data described above, we also created simulated filtered, unphased sequence data with genotype errors. We removed variants with more than two alleles or with minor allele frequency <0.5%, added genotype error at a rate of 0.1%, and removed information about genotype phase. The minor allele frequency filter of 0.5% was chosen to be higher than the threshold used to determine true IBD (0.25%). This allows the variants with frequency in the 0.25–0.5% range to be used to assess accuracy of the detected IBD segments. Error rates in sequence data vary considerably, depending on the depth of sequence coverage. We chose to use a per-variant error rate twice that of the simulated SNP data. The number of genotype errors occurring in an IBD segment is also increased in sequence data since there are more variants and thus more opportunities for error. We analyzed 500 individuals of the simulated unphased filtered sequence data, using Refined IBD with a minimum segment length of 0.2 cM and LOD scores of 4 and 5.

Wellcome Trust Case Control Consortium 2 data

The Wellcome Trust Case Control Consortium phase 2 controls consist of 5200 individuals from the United Kingdom genotyped on a custom Illumina array (Barrett et al. 2009). After quality control, including removal of SNPs not in Hardy–Weinberg equilibrium (P < 10⁻⁵) and SNPs with minor allele frequency <1%, 885,127 autosomal SNPs remained for analysis. We ran Refined IBD with a minimum genetic length threshold of 0.5 cM. We used windows of 10,000 markers (window = 10,000) with 1000-marker overlap between adjacent windows (overlap = 1000). Other parameters were left at their default values. Genetic lengths of detected IBD segments in these data and in the Finnish data (described below) were determined using the estimated genetic distances provided by the International Haplotype Map Consortium (Frazer et al. 2007).

Simulated sequence data

We generated simulated sequence data on 500 individuals. For the Refined IBD analysis, we excluded all variants with minor allele frequency ≤0.5% and added genotype error at a rate of 0.1%, which is twice that of the SNP data. Figure 6 compares IBD detection results for the simulated sequence data with those for the simulated SNP data. In the sequence data we found we needed to use a higher LOD score threshold than in SNP data to control false IBD segment discovery to a similar level. However, it should be noted that the overall level of reported segments is much higher in the sequence data, so the ratio of false to true discoveries is still well controlled with a LOD score of 3. For the same false positive IBD detection level, we detect ~50% more true IBD in the sequence data (Figure 6A). Figure 6B shows that the increase in power is due to increased ability to detect segments <1 cM. On the other hand, some of the IBD in segments >2 cM is being missed. One reason for the missed parts of long segments is the higher level of genotype error, due to both a higher error rate per variant and a higher density of variants. The phasing of the sequence data may also have a higher number of switch errors per centimorgan because low-frequency variants are generally more difficult to phase than high-frequency variants. Figure 6C shows that the accuracy of the IBD segments detected from the simulated sequence data is very high, even for segments as short as 0.25 cM.

Open in a separate window

Figure 6

Identity-by-descent detection accuracy in sequence data. Simulated sequence data on 500 individuals were used. Results from SNP data, reproduced from Figure 2, are shown for comparison. See Figure 2 for description of the axis labels.

Wellcome Trust Case Control Consortium 2 data

We analyzed ~900,000 autosomal SNPs genotyped on 5200 individuals from the United Kingdom. The average amount of IBD detected on the autosomes was 14.4 cM per pair of individuals. This equates to an IBD detection rate of 0.0041. Only 224 pairs (0.0017% of pairs) had no detected IBD. Figure 5B shows the distribution of detected IBD lengths, while Figure 7A shows the distribution of the amount of detected IBD per pair of individuals.

Open in a separate window

Figure 7

Histogram of sum of lengths of detected IBD shared by pairs of individuals. (A) In the Wellcome Trust Case Control Consortium 2 United Kingdom data. (B) In the Northern Finland Birth Cohort data. A LOD score threshold of 3 was used in both cases.

Northern Finland Birth Cohort data

We analyzed ~300,000 autosomal SNPs genotyped on 4899 individuals from Northern Finland after excluding close relatives. The average amount of IBD detected on the autosomes was 51.5 cM per pair of individuals. This equates to an IBD detection rate of 0.015, which is 3.6 times as high as that in the United Kingdom data, even though the United Kingdom data have much higher SNP density and thus much better power to detect small IBD segments. Only 675 pairs (0.0056% of pairs) had no detected IBD. Figure 5C shows the distribution of detected IBD lengths, while Figure 7B shows the distribution of the amount of detected IBD per pair of individuals.

The Northern Finland Birth Cohort (NFBC1966) Study is conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with the Broad Institute, University of California at Los Angeles (UCLA), University of Oulu, and the National Institute for Health and Welfare in Finland. This article does not necessarily reflect the opinions or views of the NFBC1966 Study Investigators, Broad Institute, UCLA, University of Oulu, the National Institute for Health and Welfare in Finland, and the NHLBI. This study makes use of data generated by the Wellcome Trust Case Control Consortium. A full list of the investigators who contributed to the generation of the data is available from www.wtccc.org.uk. Funding for the project was provided by the Wellcome Trust under awards 076113 and 085475. This study was supported by research grants HG004960, HG005701, GM099568, and GM075091 from the National Institutes of Health.