Pairs of individuals within a population can share one or more long segments of their genomes identical by descent due to inheritance from common ancestors. Identity by descent (IBD) segments are used in many applications, including estimation of kinship, 1-3 recent demography, 4-8 mutation rates, [9] [10] [11] [12] and recombination rates 13 and detection of recent selection. [14] [15] [16] The three main types of test for recent positive selection are based on population differentiation, admixture proportions, and haplotype structure. The first type looks for variants that differ markedly in frequency between populations. 17, 18 The second type looks for regions in which the sample ancestry proportions in admixed individuals differ from those elsewhere in the genome. 19, 20 One subtype of this test involves archaic admixture, such as introgression from Neanderthals into modern humans, and searches for regions in which the frequency of the archaic haplotype in a modern population is unusually high. 21 The third type looks for high-frequency haplotypes that are unusually long. 22, 23 IBD-based selection scans fall into this category. 14 IBD scans look for genomic regions that have a significantly higher than average number of IBD segments. If the genome were completely neutral, and there are no biases in detecting IBD segments or estimating their centi-Morgan (cM) lengths, the expected number of IBD segments exceeding some cM length threshold would be constant across the genome. In contrast, if certain haplotypes in a genomic region have a selective advantage, the effective size of the population is reduced in that region, which leads to a higher than expected number of IBD segments.
IBD-based tests can also detect the effects of negative selection and balancing selection, since any type of selection will tend to decrease the effective population size within the genomic region.
An IBD segment for a pair of haplotypes is a segment of DNA inherited from a single common ancestor, with no crossovers occurring within the segment in the lineages of the two haplotypes since the common ancestor. 4, 6 Within a shared IBD segment, sequence identity can be disrupted by mutation and gene conversion. In addition, genotype error can cause two haplotypes to appear to be discordant at a position. At such positions, two ''identical by descent'' haplotypes are in fact not identical. This non-identity needs to be considered when detecting IBD segments.
In the human genome, de novo single-nucleotide mutations occur at an average rate of around 1.3 3 10 À8 per base pair per meiosis, 12 which is similar to the average rate of crossing over per base pair per meiosis. Thus, regardless of the number of generations since the most recent common ancestor, an average of approximately one mutation is expected in the lineage of an IBD segment.
For a pair of haplotypes drawn at random from an outbred population, most of the genome is comprised of very short segments of IBD, with a very large number of generations to the most recent common ancestor. Since each short segment contains an average of approximately one discordance caused by mutation in addition to discordances caused by gene conversion, a series of closely spaced discordances is a clear indication that the genomic interval containing the discordances is comprised of a sequence of short IBD segments. In contrast, when one observes a long segment without discordances, it is usually (depending on the population's demographic history) highly probable that this segment is primarily comprised of a single long IBD segment resulting from recent common ancestry.
There are three primary paradigms for IBD segment detection. The first paradigm considers a pair of haplotypes to be either ''IBD'' or ''not IBD'' at each position in the genome. A hidden Markov model, with pre-determined IBD proportion and rates of transition between the IBD and non-IBD states, may be used to obtain posterior probabilities of IBD and non-IBD at each position. [24] [25] [26] [27] [28] [29] This paradigm developed out of the analysis of pedigree data and is very natural in that setting. 30 However, for population data with unknown relationships, the dichotomy into IBD and non-IBD is artificial and ignores the fact that each pair of haplotypes has a common ancestor at each position in the genome, although that ancestor may have lived a long time ago.
The second paradigm considers the length of shared segments. A segment is identical by descent if it is inherited from a common ancestor and exceeds a length threshold. In practice, if identity by state (IBS) sharing extends beyond some threshold, the segment is reported as identical by descent. [31] [32] [33] This paradigm recognizes the potential existence of IBD segments that are shorter than the threshold, but does not try to find them. The threshold is typically chosen to be a length above which the accuracy of the reported segments is high. 31, 32 The third paradigm considers two haplotypes to be identical by descent if their time to most common ancestor (TMRCA) is less than some specified number of generations. 16, 34 In this work we take a different perspective. We recognize that a pair of haplotypes is, strictly speaking, identical by descent at every point in the genome. However, for any given point in the genome, the endpoints of the IBD segment containing that point are unknown ( Figure 1A ). It is possible that one or more discordances at the end of the segment are actually contained within the long IBD segment ( Figure 1B) . It is also possible that IBD ends before IBS ends, so that the end of the IBS segment is not part of the long IBD segment, but instead contains one or more neighboring short IBD segments ( Figure 1C ). In some cases, two or more long IBD segments in a region can be mistaken for a single long IBD segment ( Figure 1D ). 35 Our approach quantifies this uncertainty. In Results, we show that our quantification is well calibrated, and we apply our method to perform an IBD-based selection scan in the UK Biobank.
The input data for our method are phased genotypes and candidate IBD segments. Highly accurate phased genotypes can be obtained from statistical phasing in large cohorts of accurately genotyped individuals. 36 The candidate IBD segments may be obtained using a length-based IBD detection method such as hap-ibd. 33 Our method estimates the posterior probability distributions of the endpoints of the candidate IBD segments and outputs quantiles and samples from these posterior distributions.
Before estimating segment endpoints, we apply a minor allele frequency (MAF) filter to the phased genotypes that excludes variants with frequency less than 0.1%. In Results we show that these rare variants are not modeled as well as the more common variants and that including these rare variants negatively impacts the accuracy of the endpoint estimates.
We model allele discordance within an IBD segment using a user-specified error rate. Analysis results are not overly sensitive to the exact choice of error rate (see Results). Discordances within a segment are assumed to occur independently except when two or more closely spaced discordances could have originated from the same gene conversion event.
We also model IBS extending beyond the end of an IBD segment. IBS segments can be comprised of multiple IBD segments. We do not try to directly model each of these IBD segments individually, but instead model the distribution of IBS segments found in the data. Short regions of IBS are modeled using the local context, because the IBS length distribution varies across the genome due to factors such as mutation rate and selection. Longer segments of IBS are modeled using chromosome-wide data because there is limited information about longer IBS segments from the local context.
We estimate the probability of the observed discordance data as a function of the IBD endpoints. We then use Bayes' rule to obtain the probability distribution of each IBD endpoint. We work from a focal position within an IBS segment ( Figure 1A ) and estimate the probability distributions for the positions of the left and right endpoints of the IBD segment that covers the focal position.
We wish to estimate the endpoints of the IBD segment covering a position x 0 for a given pair of haplotypes, H 1 and H 2 . All positions are measured in terms of genetic distance in Morgans, and haplotype phase is assumed to be known. In this description we are only concerned with the estimation of the right endpoint of the IBD segment covering x 0 . The estimation of the left endpoint is similar. Index the markers to the right of x 0 by 1; 2; 3; .; M À 1, where M À 1 is the last marker on the chromosome. Let the positions (in Morgans) of these markers be x 1 ; x 2 ; x 3 ;.;x MÀ1 . In addition, we add a nominal position x M ¼ x MÀ1 þ 1, which is located 1 Morgan beyond the last marker on the chromosome. This additional position is used to model IBD that extends beyond the end of the chromosome. We define haplotypes H 1 and H 2 to have discordant alleles at x M . The IBS data, D is the observed IBS status (identical or discordant) for the alleles on haplotypes H 1 and H 2 at the markers to the right of x 0 . Let D[a,b] denote the IBS status at markers with indices a%i%b.
Let ε be the average proportion of discordant markers within IBD segments (the ''error rate'' mentioned above). We approximate ð1 ÀεÞ with 1 and thus omit terms of ð1 ÀεÞ in our calculations.
We model the IBS data, D, from the focal point x 0 rightward as being generated by two processes. The first, up to the IBD endpoint R, requires that alleles should be identical except at a small number of discordances due to mutation, gene conversion, or genotype error. If discordances in the IBD segment are independent and the right endpoint is in the interval ðx i ; x iþ1 Þ, the probability of the data in the part of the IBD segment to the right of the focal point (i.e., PðD½1; iÞ) is e ni , where n i is the number of discordances between the first marker after the focal point and the i th marker (inclusive) and factors of ð1 ÀεÞ are approximated by 1.
We use Bayes rule to obtain the posterior distribution of R, the position of the right endpoint ðR > x 0 Þ. For each inter-marker interval ðx i ; x iþ1 Þ, i ¼ 0; 1; .; M À 1, the probability that the right endpoint is contained in the open interval ðx i ; x iþ1 Þ satisfies:
where ''f'' denotes proportionality. Normalizing the probabilities in Equation 1 to sum to 1 over the i gives the posterior probability that the endpoint occurs in each interval ðx i ; x iþ1 Þ. The probability PðR˛ðx i ; x iþ1 ÞÞ is a prior probability for the position of the right endpoint, given the position L 0 of the left endpoint (see Iterative Updating of Endpoints and Focal Point below). We model the population as having constant effective size 10,000 to obtain these probabilities. Details are given in Appendix A.
The remaining component in Equation 1 is the probability PðD½i þ1; M jR˛ðx i ; x iþ1 ÞÞ, which is the probability of the IBS data to the right of the right IBD endpoint. We model this by considering the points of discordance as being the points of renewal in a renewal process. That is, we obtain probabilities of the length of each segment that contains all of the non-discordant positions up to and including the next discordance, with each such segment being treated as being independent. The probabilities of each of these segments is obtained empirically from the observed data. Details are given in Appendix B.
Appendix C describes how to obtain the posterior cumulative distribution function for the endpoint from the interval probabilities given in Equation 1.
In the preceding section, we assumed that the left endpoint L 0 was known, but it also needs to be estimated, and we do this estimation iteratively. We start by using the left endpoint of the input candidate IBD segment as the value of L 0 . After estimating the posterior distribution of the right endpoint, we use this distribution to obtain a new ''right endpoint'' R 0 that is set equal to the 5 th percentile of this distribution. Percentiles are referenced by distance from the focal point x 0 . Thus, small percentiles are located closer to the focal point than larger percentiles. This choice of percentile is conservative (it reduces the estimated length of the IBD segment for the purpose of these calculations) and thus is likely to speed the convergence of the iterative approach. After estimating the right endpoint distribution, we estimate the left endpoint using the newly calculated value of R 0 and obtain a new value of L 0 as the 5 th percentile of the left endpoint distribution. Then if L 0 À x 0 is significantly altered (>10% change in length) from the previous value, we use the new value of L 0 to re-estimate the right endpoint distribution and obtain a new value of R 0 . If R 0 is significantly altered from the previous value, we use the new value of R 0 to re-estimate the left endpoint, and so on. Whenever we change the value of L 0 or R 0 , we update the focal point x 0 which is located half-way between L 0 and R 0 in base coordinates. We perform a maximum of 10 updates of each endpoint. In order to prevent the focal point from moving outside the input candidate IBD segment, we constrain L 0 and R 0 to stay within the input candidate IBD segment.
After running the endpoint estimation algorithm, our ibd-ends software estimates the error parameter ε by measuring the rate of discordant alleles within inferred IBD segments. The procedure is as follows. For each IBD segment that has been analyzed with the endpoint estimation algorithm, take the interval bounded by the posterior 5 th percentile of the left and right endpoint distributions. We use the 5 th percentiles so that that we are unlikely to include alleles beyond the end of the IBD segment. If the length of this interval is <2 cM, ignore the segment. Within the genomic region bounded by these endpoints, examine the alleles on the two IBD haplotypes. Count the number of mismatches and the total number of positions examined. Across all segments, report the total number of mismatches, divided by the sum of the number of positions examined in each segment. If the estimated error rate does not differ significantly from the error rate used in the analysis (e.g., less than a 3-fold difference), it is not necessary to re-run the analysis with the new value (see Results). For a large study, a pilot analysis on a small chromosome can be used to determine the error rate that should be used in the full analysis.
Gene conversions copy material from one haplotype to the other during meiosis and can thus result in discordant alleles between IBD haplotypes. The typical length of a gene conversion tract is around 300 base pairs. 37 Changes will occur only at positions at which the individual in whom the gene conversion occurred was heterozygous. Thus, many gene conversions have no effect on allele discordance, but some gene conversions can result in more than one allele discordance occurring in proximity. Since these discordances are not independent events, we do not include an error term ε for each one, since that would be overly harsh and tend to result in premature truncation of the IBD segment.
Instead, when more than one discordance occurs within 1 kb, we apply the error rate ε for the first discordance, and a less severe gene conversion error rate of ε 0 for each successive discordance within 1 kb of the first discordance (by default, e ¼ 0:0005 and e 0 ¼ 0:1).
Our software, ibd-ends, requires the input of candidate segments for which endpoints will be evaluated. In this work, we use hapibd 33 to find the candidate segments. For many applications, one wishes to assess endpoint uncertainty for all IBD segments that exceed a length threshold. In that case, the key consideration is to avoid false negatives when detecting candidate IBD segments. If a potential IBD segment is not included in the input data to ibdends, it will not be included in the results. False positives (candidates for which the true IBD segment is actually shorter than the threshold) are less serious-they increase compute time but will be shown to be unlikely to be true long IBD segments when the segment endpoints are estimated. Thus, one should try to cast a wide net when identifying candidate IBD segments.
When analyzing sequence data, the high density of variants and the presence of genotype error can cause a high rate of discordances between IBD haplotypes. The hap-ibd method permits some discordances in a segment. It does so by finding seed IBS segments that exceed a certain length, and then extending these segments if there is another IBS segment that exceeds a minimum extension length and that is separated from the seed segment by a short non-IBS gap. One way to apply hap-ibd to sequence data is to reduce the seed and extension lengths, which effectively increases the permitted density of discordances, but this can significantly increase computation time. Here we take a different approach. We reduce the marker density of the sequence data for the hap-ibd analysis (but not for the ibd-ends analysis). We choose a minimum MAF for hap-ibd that reduces the marker density to approximately that of a 600k SNP array, and we apply this threshold using hap-ibd's ''min-mac'' parameter. We retain the highest MAF variants because these are the most informative for detecting the candidate segments. This approach greatly reduces the density of variants, and thus reduces the number of IBD segments that would otherwise go undetected due to genotype error.
Except as otherwise noted, genotype data were phased using Beagle 5.1, 38 input IBD segments for ibd-ends were obtained using hap-ibd, 33 and default parameters were used for all programs.
We generated three sets of simulated data to investigate three conditions under which estimation of endpoints could be challenging: gaps in marker coverage, non-constant population size, and heterogeneous samples.
We added genotype error to the simulated data for each marker at a rate equal to the minimum of 0.02% and one-half the MAF for the marker. This error rate produces a discordance rate of 0.04% in markers with MAF > 0.04%, which matches the discordance rate seen in the TOPMed data 39 and is six times higher than the 0.0067% discordance rate seen in the UK Biobank data. 40 We also wanted to confirm that the method produces accurate results with higher error rates, so we added error to one dataset at a rate equal to the minimum of 0.1% and one-half the MAF for the marker.
True IBD segments of length Rc cM were obtained from the ancestral recombination graph. For a pair of haplotypes, we note their most recent common ancestor (MRCA) at grid-point positions located every c=10 cM along the simulated chromosome. If the same MRCA is found at R7 grid-point positions, with gaps comprising another MRCA (potentially due to gene conversion) occurring only at isolated (i.e., not adjacent) grid-point positions, we consider the pair to be identical by descent across the region. We then check all MRCAs at positions located near the first and last of the grid-points having the common MRCA to determine the exact endpoint of the IBD segment. We discard any such segments with length <c cM.
We generated 60 Mb of data for 2,000 individuals from a population with constant size of 10,000 diploid individuals. The recombination rate is 1 3 10 À8 per base pair per meiosis. During the most recent 5,000 generations, gene conversion initiations occurred at a rate of 2 3 10 À8 and gene conversion tracts had a mean length of 300 base pairs. We used SLiM v.3.3 to simulate the past 5,000 generations 41 and msprime to add mutations and simulate the more distant past. 42, 43 The mutation rate varied along the simulated region, with a new mutation rate each 100 kb that was uniformly distributed between 0 and 3 3 10 À8 per bp per meiosis. In addition, we made a 3 Mb gap by removing genetic markers between positions 20 and 23 Mb, to represent a centromeric region. We added genotype error at a rate of 0.02% as described above. We used the simulated ancestral recombination graph to determine the true endpoints of IBD segments of length >1 cM for all pairs of individuals within a subset of 500 individuals so that we could evaluate the accuracy of the inferred IBD segment endpoints.
We used the true (simulated) haplotype phase, including any alleles changed by the addition of genotype error, for all analyses of these data. When detecting candidate IBD segments with hap-ibd, we used a minor allele count threshold of 1,700 (minor allele frequency of 0.425; see Analysis Pipeline), resulting in 10,760 markers after excluding markers in the 3 Mb gap region, which corresponds to a mean density of one marker per 5.4 kb in the remaining 57 Mb. All markers with MAF > 0.1% (241,010 markers) were used in the ibd-ends analysis.
We generated 60 Mb of data for 50,000 individuals from a UK-like population. These simulated data have been described previously. 33 The demographic model has a population size of 24,000 in the distant past, a reduction to 3,000 occurring 5,000 generations ago, growth at rate 1.4% per generation starting 300 generations ago, and growth at rate 25% beginning 10 generations ago. The mutation rate is 1.3 3 10 À8 per base pair per meiosis, while the recombination rate is 1 3 10 À8 per base pair per meiosis. During the most recent 5,000 generations, gene conversion initiations occurred at a rate of 2 3 10 À8 and gene conversion tracts had a mean length of 300 base pairs. We used SLiM v.3.3 to simulate the most recent 5,000 generations 41 and msprime to add mutations and simulate the more distant past. 42, 43 We generated two copies of the data: one with 0.02% added genotype error and one with 0.1% added genotype error as described above. We also created a SNP-array version of the data with 0.02% genotype error and 10,000 randomly selected markers with minor allele frequency >5% (1 marker per 6 kb on average, corresponding to approximately 500k markers genome-wide). We used the simulated ancestral recombination graph to determine the true endpoints of IBD segments of length >1 cM for all pairs of individuals within a subset of 1,000 individuals so that we could evaluate the accuracy of the inferred IBD segment endpoints.
When applying hap-ibd to the UK-like sequence data to provide candidate IBD segments for ibd-ends, we applied a minor allele frequency threshold of 0.45 (see Analysis Pipeline), resulting in 11,524 markers across the 60 Mb with a mean density of one marker per 5.2 kb. All markers with MAF > 0.1% (198,566 markers) were used in the ibd-ends analysis of these candidate IBD segments.
We also ran hap-ibd on the UK-like sequence data with the hapibd parameters that are suggested for sequence data (min-seed ¼ 1.0, min-extend ¼ 0.2, and a minor allele frequency filter of 10%) 33 for comparison with the ibd-ends analysis.
We simulated 10 Mb of data for 500 individuals of African-like ancestry and 500 individuals of European-like ancestry. The demographic history is the two-population model of Tennessen et al., 44, 45 implemented in stdpopsim. 46 The combined sample represents an ancestrally heterogeneous population, which violates the assumption that all pairs of individuals in the sample have the same distribution of IBS segment lengths. The recombination rate and mutation rate are both 1 3 10 À8 per base pair per meiosis. We did not include gene conversion in the simulation. We simulated the data with msprime, 42 and we added genotype error at a rate of 0.02% as described above. We used the simulated ancestral recombination graph to determine the true endpoints of IBD segments of length >0.5 cM for all pairs of individuals so that we could evaluate the accuracy of the inferred IBD segment endpoints.
When applying hap-ibd to these data, we applied a MAF threshold of 0.35 (see Analysis Pipeline), resulting in 2,215 markers across the 10 Mb with a mean density of one marker per 4.5 kb. Although we used default settings for hap-ibd in other analyses (except as noted), in this analysis we set the hap-ibd minseed and min-output parameters to 1 cM since there are very few IBD segments of length >2 cM in these data. All markers with MAF > 0.1% were used in the ibd-ends analyses (48,074 markers).
We phased QC-filtered UK Biobank data (487,373 individuals) using Beagle 5.1, 38 and then used hap-ibd with default settings to find candidate IBD segments among 408,883 white British individuals identified by the UK Biobank. 40 We ran ibd-ends with default settings on the candidate IBD segments from the white British individuals to estimate the uncertainty in the endpoints of these IBD segments. We used Bherer et al.'s European genetic map which is based on family data from Iceland and other European populations. 47 Variants located outside the bounds of the map are excluded from the ibd-ends analyses because extrapolated cM positions for markers outside the map can differ significantly from their true cM positions, leading to substantial under-or over-estimation of IBD segment lengths.
Simulated Data with Variable Marker Density Figure 2 shows the results of ibd-ends analysis on the simulated sequence data with constant population size, variable marker density, and a 3 Mb gap in marker coverage (see Material and Methods for details). Even with the large gap and uneven marker density, the endpoints uncertainty is well calibrated (Figure 2A) , and coverage of IBD segments across the simulated region is even except for dips in coverage across the gap and at the chromosome ends ( Figure 2B ). 61% of sampled endpoints are located within 5 kb of the true endpoint, as are 72% of posterior median endpoints.
The ibd-ends analyses of these data used the default error rate of 0.0005. The error rate estimated by ibd-ends was 0.00039. Figure 3 shows calibration and distribution of uncertainty for the simulated UK-like data with a 0.02% genotype error rate. The estimated uncertainty is well calibrated (upper The American Journal of Human Genetics 107, 1-16, November 5, 2020 5 row of Figure 3 ), even when using inferred haplotype phase and when using data thinned to represent a 500k SNP array. Calibration is a little better with the inferredphase SNP array data rather than with the inferred-phase sequence data, which may be because phasing of common variants is more accurate than that of rare variants. As expected, average endpoint uncertainty, as measured by the difference between the endpoint sampled from the uncertainty distribution and the median endpoint, is much higher when analyzing SNP array data rather than full sequence data (lower row of Figure 3 , right versus left and middle columns). Results are also well calibrated with smaller sample sizes (200 or 1,000 individuals; Figure S1 ). We performed further analyses to evaluate endpoint estimation accuracy with a higher genotype error rate ( Figure S2 ) and with different MAF thresholds ( Figure S3 ). The results are well calibrated with the higher genotype error rate (0.1%) when using true haplotype phase. When using inferred haplotype phase, the higher genotype error rate results in phase errors that reduce accuracy. The results are not particularly sensitive to the choice of MAF, but some miscalibration is observed when very rare variants are included (0.01% MAF), and uncertainty increases when a high MAF threshold is used.
In a recent analysis, hap-ibd and GERMLINE gave the highest accuracy among competing IBD segment detection methods on sequence data, 33 so we compared the precision of estimated IBD segment endpoints between hapibd and ibd-ends for data with true haplotype phase ( Figure S4 ). Using recommended parameters for sequence data for hap-ibd (min-output ¼ 2.0, min-seed ¼ 1.0, minextend ¼ 0.2, and a minimum minor allele frequency of 10%), 33 22% of estimated endpoints were more than 50 kb from the true value when the genotype error rate was 0.02% and 43% were more than 50 kb from the true value when the genotype error rate was 0.1%; the corresponding percentages with ibd-ends were 10% and 12%. Further, hap-ibd missed many segments when the error rate was high. Only 58% of the IBD segments with true simulated length >3 cM were found by hap-ibd with the recommended parameters when the genotype error rate was 0.1%. Increasing the minor allele frequency to a high value, as we did when running hap-ibd to provide candidate segments to ibd-ends, reduces loss of IBD segments. With The data comprise 50,000 individuals simulated from a UK-like demographic history (see Material and Methods), with a genotype error rate of 0.02%. True IBD segment endpoints were determined for 1,000 individuals, and these individuals were used to generate the results in this figure. The top row shows quantile-quantile plots that assess the calibration of the estimated endpoint uncertainty. The y ¼ x line is shown for comparison. The actual quantile (y axis) corresponding to a given nominal quantile (x axis) is the proportion of segments for which the reported nominal quantile of the right endpoint is greater than the true right endpoint. The bottom row shows histograms of the right endpoint sampled from the estimated posterior distribution minus the posterior median right endpoint. The histograms represent the distribution of uncertainty, averaged over segments. Histogram bin widths are 5 kb. Results for the left endpoints are similar but are not shown. The left column is for analysis using the true haplotype phase. The middle column is for analysis using haplotype phase inferred using Beagle 5.1. The right column is for data thinned to match a SNP array with 500,000 markers genome-wide (10,000 markers in the simulated 60 Mb interval), and with haplotype phase inferred using Beagle 5.1.
hap-ibd with a high minor allele frequency threshold (45%), 88% of the segments with true length >3 cM were found when the genotype error rate was 0.1%. However, using hap-ibd with a high minor allele frequency filter without estimating the segment endpoints with ibd-ends significantly reduces precision ( Figure S4 ).
We found that using an analysis error rate that is up to three times higher or lower than the estimated error rate gave well-calibrated results ( Figure S5 ). When using an analysis error rate that is outside this range, the estimated error rates produced by ibd-ends are within the range of error rates that will provide good results in a subsequent analysis. For example, analysis of 1,000 UK-like simulated samples using an analysis error rate of 5 3 10 À5 (ten times lower than optimal) led to an estimated error rate of 4.2 3 10 À4 , while using an analysis error rate of 5 3 10 À3 (ten times higher than optimal) led to an estimated error rate of 6.2 3 10 À4 , with both estimated error rates being close to the optimal rate for these data of around 5 3 10 À4 . If the estimated error rate differs from the analysis error rate by more than 3-fold, we recommend repeating the analysis using the estimated rate. Pilot results on a small chromosome can be used to determine whether the analysis error rate needs to be changed from the default value. When we initially ran the analysis of the SNP-array data using the default error rate of 0.0005, we obtained an estimated error rate of 0.00016, which is 4-fold lower, so we re-ran the analysis with this new rate, and the results shown in Figure 3 reflect this second analysis.
Compute times for ibd-ends analyses with default settings were 1.3 h for the full UK-like data with 50,000 individuals (17 million IBD segments) and 0.5 h for 1,000 individuals (7,000 IBD segments) using a 24-core compute node with 24 Intel Xeon Silver 4214 2.2 GHz processors and 382 GB of memory.
In the heterogeneous simulation, half of the simulated individuals are from a population with an African demo-graphic history, while the other half are from a population with a European demographic history (see Material and Methods). Analyzing these data together violates the assumption that all pairs of individuals have the same distribution of IBS segment lengths, but the results are not excessively mis-calibrated ( Figure S6A ). For example, 18% of the true endpoints are closer to the center of the IBD segment than the nominal 10 th percentile. When analyzing the African individuals separately ( Figure S6B ) or the European individuals separately ( Figure S6C) , the results are well calibrated.
We determined the endpoint uncertainty distributions of 77.7 billion candidate autosomal IBD segments of length 2 cM or larger that hap-ibd detected in the 408,883 white British UK Biobank individuals. In the downstream analyses, we use the posterior medians to define the IBD segment endpoints. Analysis was parallelized by chromosome. Total wall clock computing time across all chromosomes was 7.5 h for hap-ibd and 144 h for ibd-ends using a 24-core Intel Xeon Silver 4214 2.2 GHz compute node. The estimated error rates from each chromosome varied from 0.00027 to 0.00034 when using the default analysis error rate of 0.0005. There were 54.7 billion ibd-ends IBD segments with length >2 cM. Every 10 kb along each chromosome, we computed the number of IBD segments with length >2 cM covering the position (Figure 4) . The IBD rate is the number of IBD segments covering a position divided by the number of haplotype pairs. Each individual contributes two haplotypes, and all haplotype pairs are considered except those pairs within the same individual. A high rate of IBD at a position is a signal of possible recent strong natural selection. [14] [15] [16] The median IBD rate is 0.0132%. The standard deviation of IBD rate varies by chromosome from 0.0020 to 0.019; the chromosome with the highest standard deviation is chromosome 2 because of the large spike in IBD rate around LCT (MIM: 603202) (see Figure 4) , and the median is s med ¼ 0:0027. There are 12 regions with an IBD rate higher than 0.025% (Table 1) , which is more than 43 s med higher than the median IBD rate. Eleven of the twelve regions are known to be undergoing significant levels of selection, indicating the success of this approach in finding real signals of selection. . In all cases, the genes discussed below are completely within the region for which the IBD rate is within 80% of the corresponding peak rate (column 3 of Table 1 ). The highest selection signal, with an IBD rate of 0.17%, comes from a chromosome 2 region containing LCT which has a variant selected for lactose tolerance in Europeans. 48 The selected variant is thought to have arisen, or at least begun to increase in frequency, around the time of the advent of cattle farming in Europe, around 7,500 years ago, 61 and selection has been so strong that the selected variant allele frequency is now around 75% in individuals of British descent. 62 Since IBD is affected by recent selection, it is not surprising that this signal is most prominent.
In contrast, the immunity-related HLA genes in the major histocompatibility complex (MHC) on chromosome 6 have been under selection over a much longer time period, 53, 63 and this region has a much lower peak IBD rate (0.030%) than the LCT region. Various sub-regions within the MHC appear to have been subject to adaptive introgression from Neanderthals and Denisovans, but it is difficult to be certain because long-term balancing selection across the region can produce signals that look like adaptive introgression. 64 The second-highest signal, with an IBD rate of 0.054%, comes from a chromosome 15q13 region containing TRPM1, a pigmentation gene that has been shown to have been subject to selection in non-Africans. 57, 65 The high IBD rate region on chromosome 12q24 (0.032% IBD rate) encompasses the OAS locus (OAS1-3 Regions in which the maximum IBD percentage is at least 0.025% are shown. Positions are in GRCh37 coordinates. a Region in which the IBD rate is at least 80% of the value at the peak. b See the main text for further notes and references. Genes listed within this column are contained with the region given in the third column. [MIM: 164350, 603350, 603351]) which is involved in immunity. 66 This locus has a Neanderthal haplotype present at high frequency in non-Africans that has been subject to positive selection. 56 The high IBD rate region on chromosome 17q21 (0.027% IBD rate) encompasses the 17q21.31 MAPT inversion, for which the H2 form has undergone positive selection in Europe and is associated with increased fertility and higher recombination rates in females in Iceland. 58 The high IBD rate region on chromosome 4p14 (0.027% IBD rate) contains several toll-like receptor genes (TLR1, TLR6, and TLR10) that are involved in immunity, and this region has experienced adaptive introgression from Neanderthals. 55 As well as the adaptive introgression, this region shows other signals of selection, including geographic differentiation within the UK 51 and signs of recent positive selection among non-Africans. 67 The high IBD rate region on chromosome 6q21 (0.027% IBD rate) contains the PRDM1 and ATG5 (MIM: 604261) genes. This pair of genes is associated with autoimmune diseases and cancer [68] [69] [70] and shows a signal of recent selection in HapMap data. 54 The high IBD rate region on chromosome 3p21 (0.026% IBD rate) is a region of adaptive introgression from Nean-derthals. 49, 55 It contains chemokine receptor genes, including CCR9 and CXCR6, that are involved in immunity. 71, 72 This region is also associated with COVID-19 disease severity. 73, 74 The high IBD rate region on chromosome 5p13 (0.026% IBD rate) contains SLC45A2 (formerly known as MATP), a pigmentation gene that has undergone recent positive selection in Europe. 52, 75 The high IBD rate region on chromosome 12q13 (0.026% IBD rate) contains the Type II Keratin genes, which code filament proteins that provide a major structural role in epithelial cells. 76 This region has experienced adaptive introgression from Neanderthals. 55 The high IBD rate region on chromosome 3q21 (0.025% IBD rate) contains ADCY5 (MIM: 600293). This gene is associated with birth weight 77 and type II diabetes 78 and is under long-term balancing selection. 50 The remaining locus on chromosome 22q11 (0.027% IBD rate) has not previously been highlighted as being under selection, to the best of our knowledge. This locus includes UBE2L3 (MIM: 603721) which is associated with multiple auto-immune diseases. 59 These associations make this locus a strong candidate for natural selection. 79 We also investigated the next six regions with highest IBD rates (Table 2 ; IBD rate between 0.021% and 0.025%). One of these regions is the epidermal differentiation complex locus on chromosome 1q21.3 (0.021% IBD rate), which is known to have undergone recent positive selection. 80 The variability of IBD rate across the genome is much lower for ibd-ends than for hap-ibd ( Figure S7 ). The standard deviation of the ibd-ends IBD rate is 0.0064%, whereas the standard deviation for hap-ibd is more than three times higher at 0.022%. Ibd-ends is a much better tool for investigating regions of potential selection than length-based methods because length-based methods can report artifactually high rates of IBD in some regions, such as regions containing large gaps in marker coverage. For example, IBD segment detection with four lengthbased IBD detection methods produce a 40-to 3,000-fold higher IBD rate at the chromosome 1 centromere in UK Biobank data compared to the background rate. 33 Misspecification of the genetic map can also lead to spurious regions of high IBD rate in IBD selection scans, because overestimation of genetic length will result in a larger number of IBD segments that pass the length threshold. To investigate the impact of map choice on the IBD selection analysis, we analyzed a subset of 50,000 white British individuals with three maps. We observed considerable differences in the regions of highest IBD rate when using different maps ( Figure S8 ). With the HapMap linkage-disequilibrium (LD)-based map, 86 there are 20 regions with IBD rate >0.05% (compared to 2 with Bherer et al.'s family-based map), 47 with 3 of these occurring at centromeres. Furthermore, the IBD rate has higher variability when using the HapMap map, with a standard deviation of 0.018% (compared to 0.0064% with Bherer Regions with maximum IBD rate between 0.021% and 0.025% are shown. Positions (in Mb) are in GRCh37 coordinates. a Region in which the IBD rate is at least 80% of the value at the peak. b Genes listed within this column are contained with the region given in the third column.
et al.'s map in this 50,000-individual subset). LD-based maps are known to be biased in regions of selection, with genetic distances underestimated in these regions, 87 which would lead to an apparent decrease in IBD rate in these regions. Thus, increases in apparent IBD rate with the LDbased HapMap map compared to Bherer et al.'s pedigreebased map must be due to other effects, such as unmodelled features of the HapMap data. The IBDrecomb map is designed to obtain uniform IBD rate across the genome. 13 Thus regions of likely selection such as LCT do not have high IBD rates when using this map. The analysis with IBDrecomb still has some regions of high IBD rate (21 regions with IBD rate > 0.05%), which are concentrated at chromosome ends (where inflation of map distances is known to occur 13 ) and at centromeric gaps. These regions are likely to be artifacts. [14] [15] [16] and two of these analyses also had very strong MHC signals. 14,16
We presented a method for calculating the posterior probability distributions of IBD segment endpoints. We showed that the method can be applied to large datasets, such as the UK Biobank SNP array data on 408,883 white British individuals and simulated sequence data on 50,000 individuals. In the UK Biobank data, we analyzed 77.7 billion candidate IBD segments and found 54.7 billion IBD segments for which the length based on posterior median endpoints is greater than 2 cM.
In addition to quantifying endpoint uncertainty, a major advantage of our method is that it handles genotype errors and other discordances within IBD segments in a principled way, in contrast to many other methods for IBD segment detection which use ad hoc approaches. Our method does not directly account for haplotype phase uncertainty, but statistical phasing of non-rare variants in large array-typed cohorts is now extremely accurate, 36 and technologies for generating highly accurate, phaseresolved sequence data are becoming available. 88 Unlike most existing methods for IBD segment estimation, ibd-ends is designed for sequence data. We found that a length-based IBD detection method could not simultaneously detect IBD segments with high power and precisely determine endpoints, but the ibd-ends method detects nearly all segments exceeding a specified length and accurately determines the endpoints, even in the presence of relatively high rates of genotype error.
IBD segments defined by the sampled endpoints or posterior quantiles such as the median can be used in down-stream analyses. This can enable the inclusion of more data in analyses. For example, when estimating genomewide mutation rates from IBD segments, it is important to be confident that one does not count mutations that actually lie outside the IBD segment. In the past, this has been achieved by trimming 0.5 cM from the putative IBD endpoints, 10,12 but trimming using a small quantile of the uncertainty distribution would be expected to result in less trim (and hence more data), while maintaining accuracy. Alternatively, methods could be developed that directly account for endpoint uncertainty without trimming.
Another analysis that could incorporate estimated IBD segment end point uncertainty is IBD-based estimation of recombination maps. Endpoints of IBD segments are points of past recombination which are used to estimate the map. Misspecification of IBD endpoints adds noise and bias to the resulting map. Thus, the higher precision of posterior median endpoints relative to endpoints from other IBD detection methods could improve estimation of recombination rate. One could also improve recombination map estimation by incorporating endpoint uncertainty into an iterative procedure that uses the current estimated recombination map to refine the estimates of IBD segment endpoints.
IBD segment endpoint uncertainty could also be used to improve IBD-based estimation of recent demographic history. For this application, it is not the actual IBD endpoints that matter, but rather the distribution of IBD segment lengths, which one could estimate by sampling endpoints from the posterior distribution. A threshold of 2 or 3 cM on IBD segment length is usually applied with current methods since lengths of shorter segments have higher relative uncertainty. [4] [5] [6] [7] [8] With sampled endpoints, it may be possible to use a lower length threshold and thus estimate demographic history further into the past.
We found that the IBD segments obtained from sampled endpoints provide excellent input data for an IBD-based selection scan. Eleven of the twelve top regions in our UK Biobank analysis are regions of known selection, and the remaining region is a good candidate for selection. Two particular features of our IBD-based selection analysis contribute to its success. The first is accurate estimation of IBD segment lengths, even in the presence of gaps in marker coverage, which eliminates many spurious signals. In contrast, length-based IBD detection methods tend to have regions with inflated IBD rates, 33 which would cause spurious signals if used in a selection scan. The second is the choice of genetic map. We found that for IBD-based selection scans, pedigree-based maps based on actual observations of recombination produce more accurate results than maps based on LD or on IBD sharing, because effects of selection and other unmodelled features in local genomic regions can bias LD-based and IBD-based maps, and this bias results in spurious signals of selection.
The modeling underlying our method assumes a homogeneous population with the same distribution of IBS segment lengths for all pairs of individuals. Analysis of a simulated combined African and European ancestry sample indicates that violation of this assumption introduces a little mis-calibration in the posterior endpoint probabilities. One solution for future work would be to adjust for ancestry of the samples, or for local ancestry in the case of admixed samples.
In this appendix we derive a formula for the prior probabilities PðR˛ðx i ; x iþ1 ÞÞ in Equation 1. Let Y ¼ R À L 0 be the length (measured in Morgans) from the left endpoint, L 0 , of the segment to the right endpoint R.
Here we assume that the left endpoint is known, but in practice we iteratively update its estimated value as described in the main text. We write FðyÞ ¼ PðY %yÞ for the prior probability distribution on Y. We model the population size as constant diploid size N, with N ¼ 10,000, which reflects the approximate average historical size of out-of-Africa populations.
Summing over possible values for G, the number of generations to the common ancestor of the IBD segment, we obtain
In these calculations, we use the fact that the IBD segment length conditional on the number of generations g to the common ancestor is exponentially distributed with rate parameter 2g, 7 and the fact that the number of generations to the most recent common ancestor in a population of constant size n is geometrically distributed. 4 The final expression for FðyÞ in the calculation is obtained by applying the formula for a geometric series.
The prior distribution for the right endpoint is:
and the probability in Equation 1 that R falls in a certain interval can be calculated as:
In Appendix C, we will require the inverse of F, which is
Appendix B: Modeling the IBS Data beyond the IBD Segment In this appendix we describe how to calculate the conditional probability PðD½i þ1; MjR˛ðx i ; x iþ1 ÞÞ in Equation 1. This is the probability of the IBS data to the right of the right IBD endpoint, including a nominal discordant posi-tion indexed by M that is 1 Morgan beyond the last marker on the chromosome. Let m i ðjÞ ðj R1Þ be the ordered indices of the discordant markers to the right of x i . The final value of this sequence is index M. We approximate the IBS process to the right of the IBD endpoint as a renewal process with a renewal every time there is a discordant marker. Then the probability of the IBS data to the right of x i given that the IBD endpoint occurred in the interval ðx i ; x iþ1 Þ is:
Each interval of IBS data in the preceding equation has the form D½a; b, where a and b are marker indices, and has the property that the alleles on the two haplotypes H 1 and H 2 are identical at all positions i satisfying a%i < b, but are discordant at position b. Our approach to estimating the probabilities of these IBS interval data is to estimate the probability p a,b that two randomly chosen haplotypes have identical alleles over the interval [a,b] (i.e., for all marker indicies i with a%i%b). Then we set PðD½a; bÞ ¼ p a;ðbÀ1Þ À p a;b : In the case where a ¼ b (due to two consecutive discordances), we define p a;ðbÀ1Þ ¼ 1, so that PðD½b; bÞ ¼ 1 À p b;b which is the probability of discordance at marker b. We estimate p a,b empirically from the data.
Our method for estimating p a,b depends on the length of the interval [a,b]. For short intervals, we use data in the interval so that the estimated probability incorporates the local genomic context, such as high or low marker density, high or low heterozygosity, and high or low levels of LD. For long intervals, we will use chromosome-wide data. As the interval length becomes longer, p a,b will tend to decrease because the probability of observing a long IBS interval for a random pair of haplotypes is small. Small probabilities are more difficult to estimate than long ones, so we need to bring in data from the rest of the genome. Furthermore, long intervals represent long IBS (alleles at all markers in the interval are identical except for the final marker), and long IBS is likely to be the result of a long IBD segment. Excluding the effects of selection, the distribution of IBD lengths should be uniform across the genome, so estimating the frequency of such long IBS segments from data across the genome is appropriate.
For short intervals, we estimate p a,b by the proportion of pairs of haplotypes that have identical by state alleles for markers in the interval [a,b]. We consider an interval [a,b] to be a short interval if the estimated b p a;b satisfies b p a;b R0:001 when it is estimated from the haplotypes in the interval. We use 10,000 randomly sampled haplotypes (or all haplotypes if the sample size is 5,000 or fewer individuals) to estimate the short interval IBS probabilities.
For longer intervals, we estimate p a,b using the global distribution of sampled one-sided IBS lengths. A one-sided IBS length is the distance in Morgans from a random starting position to the first non-IBS marker in the direction toward the center of the chromosome for a random pair of haplotypes. We estimate p a,b by the proportion of sampled one-sided IBS lengths that are greater than ðx b À x a Þ. When estimating the global IBS length, we randomly select 1,000 positions and randomly sample 2,000 one-sided IBS lengths for each position. We then exclude the sampled lengths at positions for which the IBS length are significantly longer than average. We do this by calculating the 90th percentile of the 2,000 segment lengths for each position. If the 90 th percentile at a position is more than 3 times the median 90 th percentile from all 1,000 positions, we discard the position. This filtering protects against selecting positions near gaps in marker coverage.
In this appendix we describe how to calculate the cumulative distribution function of the posterior distribution for the right endpoint of the IBD segment, as well as how to sample from this posterior distribution.
We assume that the data D are not informative about the location of the endpoint within the inter-marker interval given that the endpoint occurs within the interval, i.e., we assume that Pðx i < R % xjx i < R < x iþ1 ; DÞ ¼ Pðx i < R % x jx i < R < x iþ1 Þ for x i < x < x iþ1 , since there are no data within the interval. Write p i ¼ PðR %x i jDÞ, which can be estimated using Equation 1 with the procedures described above. As in Appendix A, we write Y ¼ R À L 0 for the length of the segment (in Morgans), and FðyÞ ¼ PðY %yÞ for the prior on IBD lengths. For x satisfying x i < x < x iþ1 : Note that p i ¼ PðR %x i jDÞ is conditional on the data, whereas Fðx i ÀL 0 Þ ¼ PðY %x i ÀL 0 Þ ¼ PðR %x i Þ is a prior probability and is not conditional on the data.
Then to find the p th quantile, we want to find x ðpÞ such that PðR < x ðpÞ DÞ ¼ p. Solving the above equation for
x ¼ x ðpÞ , we obtain
The formula for F À1 ðpÞ can be found in Appendix A.
In order to obtain a sampled value from the posterior probability of the endpoint, we first generate a realization u from the Uniform(0,1) distribution, and we then obtain x ðuÞ using the above formula with u in place of p, which is then a realization from the desired distribution. This is an example of using the inverse transform principle to sample from a distribution. 89 
The ibd-ends software is freely available under the open source Apache License 2.0 (see Web Resources).
The UK Biobank IBD rates shown in Figure 4 are available from Mendeley Data https://doi.org/10.17632/ m4nxyv6rw8.1
Supplemental Data can be found online at https://doi.org/10. 1016/j.ajhg.2020.09.010.
Research reported in this publication was supported by the National Human Genome Research Institute of the National Institutes of Health under award numbers HG005701 and HG008359. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This research has been conducted using the UK Biobank Resource under Application Number 19934.
The authors declare no competing interests.
Fðx À L 0 Þ À Fðx i À L 0 Þ Fðx iþ1 À L 0 Þ À Fðx i À L 0 Þ À p iþ1 À p i Á