{'fraction_non_alphanumeric': 0.02716198419666374, 'fraction_numerical': 0.04140144863915716, 'mean_word_length': 4.918967197141929, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identifying non-identical-by-descent rare variants in population-scale whole genome sequencing data

The site frequency spectrum in human populations is not accurately modeled by an infinite sites model, which assumes that all mutations are unique. Despite the pervasiveness of recurrent mutations, we lack computational methods to identify these events at specific sites in population sequencing data. Rare alleles that are identical-by-descent (IBD) are expected to segregate on a long, shared haplotype background that descends from a common ancestor. However, alleles introduced by recurrent mutation or by non-crossover gene conversions are identical-by-state and will have a shorter expected shared haplotype background. We hypothesized that the expected difference in shared haplotype background length can distinguish IBD and non-IBD variants in population sequencing data without pedigree information. We implemented a Bayesian hierarchical model and used Gibbs sampling to estimate the posterior probability of IBD state for rare variants, using simulations to demonstrate that our approach accurately distinguishes rare IBD and non-IBD variants. Applying our method to whole genome sequencing data from 3,621 individuals in the UK10K consortium, we found that non-IBD variants correlated with higher local mutation rates and genomic features like replication timing. Using a heuristic to categorize non-IBD variants as gene conversions or recurrent mutations, we found that potential gene conversions had expected properties such as enriched local GC content. By identifying recurrent mutations, we can better understand the spectrum of recent mutations in human populations, a source of genetic variation driving evolution and a key factor in understanding recent demographic history.


Introduction 20 21
Recurrent mutations are repeated mutational events at the same nucleotide position in multiple individuals 22 in a population. The frequency of recurrent mutations and their relevance to evolutionary genetics studies 23 have been examined since the beginning of the field of population genetics (e.g. Wright 1931;Haldane 24 1933;Wright 1937). The frequency that a recurrent mutation is observed in a sample depends on many 25 factors, including the per-base-pair mutation rate, the number of chromosomes surveyed, the effective 26 population size, as well as the demographic history of the population surveyed. Distinguishing recurrent 27 mutations from variants whose alleles are all inherited identical-by-descent (IBD) is critical to a complete 28 understanding of the human germline mutation rate, and to population genetic methods that make 29 inferences from the observed number and frequency of genetic variants in a population. 30 As the genetics community has surveyed large, rapidly growing populations with a finite genome 31 size, such as modern humans (Harpak et al. 2016), it has been observed that recurrent mutations occur at 32 appreciable frequency. For example, in the Exome Aggregation Consortium dataset of 60,706 human 33 exomes, there is a marked absence of singleton CpG transitions relative to other mutation types (Lek et al. 34 2016). This observation could be explained by the presence of recurrent mutations saturating these highly 35 mutable sites in this large sample, resulting in two or more sampled individuals segregating identical-by-36 state alleles at CpG sites. 37 One implication of this observation is that the presence of recurrent mutations in a large sample 38 may result in a suboptimal calibration of summary data typically utilized in population genetic inference, 39 like the site frequency spectrum (SFS). The SFS is the distribution of the number of observed variants at 40 allele counts 1 to n-1 in a sample of n chromosomes. Many modern population genetics methods use the 41 SFS to infer the demographic history of a sample (Gutenkunst et al. 2009;Lukic and Hey 2012;Excoffier 42 et al. 2013;Bhaskar et al. 2015;Jouganous et al. 2017). These approaches generally assume an infinite 43 sites model with no recurrent mutations, but the human site frequency spectrum is not well explained by 44 an infinite sites model (Harpak et al. 2016). Previous work has described the SFS allowing for recurrent 45 mutations, relying on observed recurrent mutations in the form of triallelic sites (Jenkins and Song 2011;46 Jenkins et al. 2014;Ragsdale et al. 2016). If recurrent mutations are not accounted for, the SFS will be 47 shifted to higher allele frequencies relative to the SFS that incorporated recurrent mutations at lower 48 frequencies. This could impact the accuracy of demographic parameter inference. In particular, methods 49 that infer the magnitude of recent population growth, which rely on rare variants, may incur bias if they 50 do not take recurrent mutations into account. Similarly, the magnitude of purifying selection may be 51 underestimated if the frequencies of rare variants are overestimated due to undetected recurrent mutation. 52 We propose to identify rare variants that fall on the extreme short end of shared IBD segment lengths 86 and attempt to categorize these likely non-IBD variants by the possibilities enumerated above. In addition 87 to the IBD segment length, additional genomic annotations can help to distinguish between these causes, 88 such as the mutation's sequence context (e.g. a CpG mutation), the local recombination rate, and local GC 89 content. 90 While previous efforts have leveraged IBD tracts to infer mutation and gene conversion rates 91 genome-wide scale. In today's era of whole genome sequencing of thousands of individuals from a 95 population, categorizing specific rare variants as likely recurrent mutations or gene conversions is now 96 uniquely possible. Here, we describe a Bayesian hierarchical model to identify non-IBD rare variants, 97 using population genetic simulations to assess its precision and accuracy. We then apply our approach to 98 sequencing data of 3,621 individuals from the UK10K dataset, and partition high-confidence non-IBD 99 rare variants as those likely to be recurrent mutations or gene conversions. Previous work to model the expected TMRCA between two IBD alleles or two random alleles in a 105 population provides a framework through which these states can be distinguished in data. Measuring the 106 accumulation of mutations on a haplotype, i.e., the mutational clock, is useful for estimating the age of 107 older, common variants; however, for our purpose here to distinguish rare recurrent and IBD variants, 108 there will be few if any linked mutations more recent than the focal variant. Therefore, the mutational 109 clock does not help us distinguish IBD and non-IBD rare alleles, and we rely solely on the recombination 110 clock for inference. 111 The theoretical distributions of the pairwise TMRCAs for IBD or non-IBD alleles for a range of 112 allele counts are plotted in Supplementary Figure 1 (Supplementary Methods). As the IBD allele 113 count increases, the mean TMRCA also increases, reflecting that higher frequency alleles tend to be 114 older; meanwhile, the TMRCA distribution between non-IBD allele pairs is unchanging because it is not 115 a function of the allele frequency. Thus, the difference between the expected TMRCA for IBD vs. non-116 IBD variants increases with decreasing allele frequency. 117 Though the TMRCA of a genetic variant is not directly observable, it can be estimated by the 118 length of the haplotype shared by carriers of the variant. The distance to the nearest recombination event 119 on either side of a genetic variant between a pair of alleles can be modeled as exponentially distributed 120 with rate proportional to the TMRCA (Palamara et al. 2012;Mathieson and McVean 2014). The expected 121 difference in the TMRCA between a rare variant segregating with IBD alleles and a variant position with 122 recurrent mutations (Supplementary Figure 1) translates into IBD variants having, on average, longer 123 pairwise distances to obligate recombination events compared to recurrent sites of the same allele 124 frequency (Supplementary Figure 2) develop an approach to identify non-IBD variants that scales to large, whole genome population 133 sequencing studies with thousands of individuals. 134 While recombination breakpoints cannot be directly observed in population sequencing data, 135 patterns of genetic variation can give us an estimate of the location of these events. Here, we are 136 interested in rare genetic variants that are difficult to accurately phase. Additionally, the signature of 137 recurrent mutation itself could introduce error into statistical phasing algorithms. Thus, our method 138 utilizes unphased diploid genotypes to estimate the recombination distances on either side of a pair of 139 alleles. With diploid genotypes for a pair of individuals each carrying a focal allele, one can measure the 140 obligate recombination distance as the physical span to the first opposite homozygote genotype between 141 the two individuals (Supplementary Figure 3). No genealogy without recombination is compatible with 142 the observed genotypes of these two sites (the focal allele and the site of the opposite homozygote 143 genotypes), and so we assume a recombination event has occurred between them (Mathieson and 144 McVean 2014). Thus, the obligate recombination distance gives an estimate of the true recombination 145 distance. 146 147

Statistical Model 148
We considered a Bayesian hierarchical model for the pairwise recombination distances from a sample of 149 variants of a given allele count, which allowed us to learn the model parameters directly from the data 150 (Figure 1). We modeled the sampled variants as a finite mixture of IBD (k=1) or non-IBD (k>1, with 151 each possible partition of alleles for a non-IBD variant given a different value of k), with mixture 152 proportions π k . For example, a non-IBD variant of allele count 4 has two possible partitions: a singleton 153 and an IBD tripleton (1:3), or two doubletons (2:2). Each variant of allele count A has = ! ! allele 154 pairs. The TMRCA (t) for an allele pair was sampled from a gamma distribution with shape α and rate β, 155 with one gamma distribution of t for IBD allele pairs and one for non-IBD allele pairs, For non-IBD allele 156 pairs, we estimated α and β from multiallelic sites, and for IBD allele pairs we fixed α and performed 157 sampling for β, over a range of possible values for α. For each variant, the possible permutations of allele 158 pair assignments of IBD or non-IBD states are denoted by j. For IBD variants, all allele pairs are IBD; for 159 non-IBD variants, the possibilities depend on the partition k. We modeled the left and right recombination 160 distances (d L , d R ) for each allele pair following an exponential distribution with rate proportional to t. We 161 used Gibbs sampling to sample from the marginal posterior density of each parameter, as we could 162 estimate these densities from the full conditional distributions. Below we outline these expressions. 163

Mixture proportions (π):
Using a multinomial likelihood for the probability of the assignments k based 164 on proportions π, we used the conjugate prior Dirichlet distribution to get a Dirichlet posterior for the 165 probabilities of π given the observed k assignments. Thus we have the likelihood function: 166 We used a Dirichlet prior for π: 167 ~ℎ ( 2 )

168
The resulting posterior probability followed a Dirichlet distribution: 169 170

TMRCA (t):
The likelihood of the pairwise recombination distance d to one side of a variant (in 171 j, the partition of non-IBD variants: 199 200

202
Application to simulated genetic data 203 To evaluate our approach, we applied it to simulated genetic data including non-IBD (recurrent) We next performed a battery of sensitivity studies, simulating population genomics features 219 known to influence patterns of genetic variation that may impact the robustness of our estimates. First, we 220 performed simulations of genomic segments including genes and deleterious mutations, and applied our 221 approach to these simulations to test the effect of selection on our approach to identify non-IBD variants 222 (Methods). We found that including background selection had little impact on the power of our approach 223 to identify recurrent mutations (Supplementary Figure 6; Figure 7); thus, the difference in recombination distance between IBD and recurrent 226 variants in these simulations is not strongly altered by the presence of weak negative selection. 227 Next, we evaluated the performance of our approach for simulated variants flanked by 10,000 228 base pair recombination hotspots, with hotspot recombination rates of 5x10 -6 , 1x10 -6 , or 5x10 -7 events per 229 base pair per generation (Methods). For variants close to a hotspot, we expected a smaller difference 230 between recurrent and IBD allele pairs' recombination distances, due to a weaker relationship between 231 recombination distance and TMRCA. When we applied our Bayesian hierarchical model to the 232 simulations with hotspots, we found that our power to distinguish IBD and recurrent variants decreased 233 non-IBD variants, and that runtc's performance at this task -a task we note that it was not explicitly 251 designed for -performs poorly compared to our Bayesian hierarchical approach (Supplementary Table  252 2). 253 254

Application of Bayesian hierarchical model to UK10K sequencing data 255
We applied our method to identify non-IBD variants in whole-genome sequencing data in 3,621 256 individuals from the UK10K project (Walter et al. 2015). Individuals from the ALSPAC and TWINSUK 257 studies used here were sequenced to average depth ~7x and passed the UK10K project quality control 258 filters. We measured the obligate recombination distance for biallelic and multiallelic single nucleotide 259 variants that passed the UK10K quality filters. Based on the decreased performance of our method with 260 increased allele count, we restricted our analysis to variants of allele count less than or equal to 5. 261 We applied our approach to a mixture of 80% biallelic and 20% multiallelic sites, in order to use 262 multiallelic sites as a positive control for non-IBD mutations. We compared the empirical cumulative 263 distribution of posterior probabilities for multiallelic and biallelic sites, and as expected we observed that 264 multiallelic sites had higher posterior probabilities of being non-IBD (Supplementary Figure 11). We 265 used these distributions to determine the threshold of posterior probabilities we called "likely non-IBD" 266 for all biallelic variants at allele counts 2-5, which we then used in downstream analyses. 267 268

Non-IBD variants correlate with local sequence context 269
To assess the accuracy of our recurrent mutation calls, we took advantage of the relationship between 270 local sequence context and mutation rate (Aggarwala and Voight 2016). Under a Poisson model of 271 mutation, sequence contexts with a higher mutation rate should have a higher probability of recurrent 272 mutation relative to other contexts (i.e., "double hits"). If non-IBD variant calls reflect recurrent 273 mutations, we would expect to see a correlation between the fraction of non-IBD variants and the 274 mutability of sequence contexts. Conversely, if our approach randomly selects a subset of sites rather than 275 true recurrent mutations, we would not expect to see a relationship between sequence context and fraction 276 of sites called recurrent. Using sequence-context estimated polymorphism probabilities calculated from 277 the UK10K dataset, we calculated an expected fraction of recurrent variants for each 5 base-pair (5-mer) 278 sequence context and allele count (Methods). 279 Across all 5-mer sequence contexts, we observed a significant correlation between expected and 280 observed fractions (e.g. Pearson's correlation = 0.81, P < 10 -100 for allele count 2; Figure 3; 281 Supplementary Table 3). The observed fraction of non-IBD called sites was higher than expected for 282 non-CpG->T contexts, and lower than expected for CpG->T contexts (Figure 3; Supplementary Table  283 4). Within CpG->T contexts, we also observed a significant correlation between expected and observed 284 fractions, though for all contexts the observed fraction of non-IBD calls was less than expected 285

Additional genomic annotations correlated with non-IBD variants 293
Next, to understand which genomic features in addition to local mutation rate are associated with non-294 IBD variants, we performed a linear regression with the posterior probability of each variant being non-295 IBD as the response variable (6,763,324 sites; with 665,340 called non-IBD). In separate regressions for 296 each allele count, we included 7-mer polymorphism probabilities, background selection, GC content, 297 replication timing, local recombination rate, distance to a recombination hotspot, germline CpG 298 methylation levels, the variant calling quality measure VQSLOD, and read depth as predictor variables. 299 We transformed the values of each annotation to z-scores, and report the odds ratios and 95% confidence 300 interval for each annotation in Figure 4 (Supplementary Table 5). All annotations were significantly 301 associated with the outcome (P<1x10 -10 ). In addition, we performed a logistic regression with IBD/non-302  Table 5). Below, we highlight the annotations included as predictors, our prior hypotheses about their 305 relationships with recurrent mutations, and the results of the regression models. 306 307 Polymorphism probability: As shown in our analysis of expected vs. observed recurrent fraction for 5-308 mer sequence contexts above, polymorphism probability was strongly positively correlated with non-IBD 309 status. As previous work has shown that a 7-mer model explains additional variation in genetic variation 310 over a 5-mer model (Aggarwala and Voight 2016), we find that a 7-mer polymorphism probability 311 calculated in UK10K in the logistic regression model was associated with our non-IBD calls. neutral variation in a region, i.e., stronger background selection. We expected that increased background 324 selection would be associated with increased recurrent mutation, as linkage to deleterious alleles would 325 result in variants being removed from the population and thus present at lower frequencies. Recurrent 326 mutations would then be more likely to be present as they effectively shift the site frequency spectrum 327 towards more rare alleles. Our results are consistent with this expectation, with an odds ratio less than one 328 for B-values. 329 Local recombination rate and distance to recombination hotspots: The results of our simulations 330 suggested that we have lower power to identify recurrent variants located near a recombination hotspot 331 (Supplementary Figures 8, 9). Indeed, we observed that both an increased local recombination rate and a 332 shorter distance to a recombination hotspot were correlated with a lower probability of a site being called 333 as recurrent. 334

Methylation levels at CpG sites: Spontaneous deamination of 5-methylcytosine at CpG sites results in a 335
substantial increase in C-to-T transition mutations. We included CpG methylation levels measured in 336 testes and ovaries in our model, expecting that CpG sites with higher methylation levels are more likely to 337 spontaneously deaminate, increasing mutation rates generally and thus increase recurrent mutation 338 probabilities. Methylation levels in testes and ovaries were correlated (Pearson's correlation coefficient = 339 0.27, P<2x10 -16 ), but we noted that increased methylation in both tissue types independently predicted an 340 Limiting our results to tracts less than 1kb with 2 or more non-IBD variants present in the same 360 individuals, we identified 42,203 variants within 18,971 putative gene conversion tracts, representing 361 6.3% of non-IBD variants (Supplementary Figure 16). We performed logistic regression with all non-362 IBD variants labeled as potential gene conversions or not as the outcome, and the genomic annotations 363 listed above as predictor variables (Figure 5, Supplementary Table 6). We additionally included the 364 posterior probability of a variant being non-IBD as a predictor variable. Compared to non-IBD variants 365 not in putative gene conversion tracts, these variants were associated with lower polymorphism 366 probability, higher variant quality score, increased posterior probability of being non-IBD, smaller 367 distance to a recombination hotspot, and lower recombination rate. We also observed a GC bias in 368 putative gene conversion variants, as measured by the fraction of variants containing an A->C/T->G or A-369 >G/T->C mutation (37% in putative gene conversions vs. 27% in all other non-IBD variants; P < 10 -100 ; 370 Fisher's exact test). 371 372

Rescaling the site frequency spectrum with recurrent mutations 373
With our set of high-confidence non-IBD variants, we rescaled the site frequency spectrum for very rare 374 variants. Taking into account the power of our approach on simulated data, we plot the original and 375 rescaled SFS for variants with allele count <5 in Figure 6A (Methods). Rescaling the site frequency 376 spectrum resulted in a 3% increase in the fraction of singleton variants, from 46.6% to 49.6%. As 377 expected, the majority of this shift is due to the relatively large fraction of CpG->T variants that were 378 called as non-IBD ( Figure 6B). For CpG->T variants alone, the fraction of singleton variants increased 379 from 43.6% to 49.9%. We note that this rescaling is incomplete, as we identified non-IBD variants at only 380 allele counts 2 to 5 (representing 38% of non-singleton variants in UK10K). 381 382

Discussion 383
We describe a novel approach designed to specifically identify non-IBD variants in whole genome 384 sequencing data by leveraging the difference in the obligate recombination distance between rare IBD and 385 non-IBD variants. Our approach uses a Bayesian hierarchical model and Gibbs sampling to jointly infer 386 the TMRCA distributions of these two scenarios and identify variants with a high posterior probability of 387 being non-IBD. In simulated data, we find that the posterior probabilities of a variant being non-IBD can 388 discriminate between IBD and recurrent mutations for variants up to allele count 5 in a population sample 389 of 3,621 individuals. 390 Our approach assumes that we do not have phase information for individuals, i.e. we do not 391 assign each variant in an individual to a maternally or paternally inherited chromosome. If we had 392 accurate phase information for rare variants, such as from long-read sequencing data, or 'hard-phase' calls 393 from paired-end sequencing libraries, we could more accurately measure recombination breakpoints. This 394 could potentially improve the accuracy of our method by eliminating the measurement error caused by 395 using the obligate recombination distance. 396 We focused on identifying non-IBD variants for allele counts of 5 or less, as the performance of 397 our method decreases with increasing allele count. Additionally, the computational burden of sampling 398 from the marginal posterior distributions increases exponentially with increasing allele count. With a 399 larger sample size, the frequency of a variant at a given allele count will decrease, while the 400 computational complexity remains the same. Thus, we expect that applying our method to even larger 401 sequencing datasets will improve its performance. 402 The negative correlation we observed between local recombination rate and the probability of a 403 site being called non-IBD suggests that our method is confounded by local recombination rate. In 404 simulated data, we also observed that we had lower power to identify recurrent mutations in close 405 proximity to recombination hotspots. We also note that we found a significant relationship between 406 sequencing quality, measured by read depth and variant quality score, and the probability of a site being 407 called non-IBD. The signature of a non-IBD variant used here could also be that of a genotyping error, as 408 genotyping errors also may occur on any random haplotype background in a population. This could be a 409 potential application of our method, as a way to identify genotyping errors in large scale sequencing 410 datasets. Our current recommendation to overcome this issue is to remove variants of low quality until 411 this relationship is not significant. However, distinguishing genotyping errors from true non-IBD variants 412 remains an important problem.  For forward genetic simulations with selection, we generated 10Mb genomic segments using a 429 recipe from the SLiM manual (Haller and Messer 2017) with the following procedure: 1) sample non-430 coding region; 2) sample exon; 3) sample intron and exon pairs in a loop with 20% probability of 431 stopping after each pair; 4) repeat steps 1-3 while chromosome length < 10Mb; 5) sample final non-432 coding region. Exonic mutations were synonymous or non-synononymous at a ratio of 1:2.31, and 10% of 433 non-synonymous mutations were neutral. Deleterious non-synonymous mutations' selection coefficients 434 were sampled from a gamma distribution with mean -0.03 and shape 0.207. Exon lengths were sampled 435 from a lognormal distribution with mean log(50) and standard deviation log(2). Non-coding regions were 436 neutral and their lengths were sampled from a uniform distribution between 100 and 5000. Intronic 437 mutations were neutral and intron lengths were sampled from a lognormal distribution with mean 438 log(100) and standard deviation log(1.5). 439

UK10K dataset 441
We applied our method to identify recurrent mutations in whole-genome sequencing data in 3,621 442 individuals from the UK10K project (Walter et al. 2015). These individuals were sequenced to average 443 depth 7x, passed the UK10K project quality control filters, and come from the ALSPAC and TWINSUK 444 studies. We measured the recombination distance for biallelic and multiallelic single nucleotide variants 445 that passed the UK10K quality filters that were present at allele count ≤ 10 in these individuals. 446 447

Measuring the obligate recombination distance 448
For simulated data, we generated diploid genotypes by randomly combining pairs of haploid genomes, 449 and calculated the recombination distances for variants within the central 2Mb of each 10Mb genomic 450 segment. In both simulated and UK10K data, we measured the obligate recombination distances for 451 variants with allele count ≤ 10. For each pair of carriers, we identified the nearest variant upstream and 452 downstream with opposite homozygote genotype, i.e. where one individual has genotype 0 and the other 453 has genotype 2 (Supplementary Figure 3). We then converted the physical distance to a genetic distance 454 using a genetic map. For UK10K, we used a 1000 Genomes Project CEU genetic map (The 1000 455 Genomes Project Consortium 2012), and for simulated data we used a uniform map with r = 1x10 -8 events 456 per site per generation, or a variable map for simulations including recombination hotspots. 457 458

Applying the Bayesian hierarchical model 459
To apply our model to simulated or UK10K recombination distances, we first generated an estimate of the 460 beta parameter for non-IBD variants from multiallelic sites. Using Gibbs sampling on non-IBD allele 461 pairs from multiallelic variants, we used a simplified version of the hierarchical model where we sampled 462 the TMRCA for each allele pair and the beta parameter in each Gibbs iteration. We repeated this 463 procedure to estimate beta for a range of alpha values from multiallelic sites' recombination distances. To 464 test if the choice of alpha affected our ability to discriminate IBD and non-IBD variants, we applied the 465 model with different non-IBD alpha/beta values to UK10K variants on chromosome 22. The posterior 466 estimates of k were highly correlated across values of alpha (alpha=20 vs. alpha=40, Supplementary 467 Table 7). When applying the full model to data, we used alpha = 10 for IBD allele pairs, with alpha = 40 468 and the corresponding value of beta inferred from multi-allelic sites (beta = 0.0859) for non-IBD allele 469 pairs. We ran 10,000 iterations of the Gibbs sampler for each run of the model, thinned the chains until 470 autocorrelation was below 0.01, and assessed convergence of the chains by comparing the thinned 471 samples from the first and second half of the chain via a Wilcoxon rank-sum test. A chain was determined 472 to have converged if the Wilcoxon test P-value was > 0.05. 473 We parallelized the application of our model by breaking down the genome into 10Mb segments, with thousands of variants in each segment, we observed no effect, but for larger allele counts we did see 478 an effect of applying the model to small numbers of variants. Thus, for allele counts >5, we grouped 479 segments together until at least 1000 variants were included in each run of the model. For all ROC curves from simulated data, we calculated the area under the curve as: 489 where r and i represent recurrent and IBD variants, and q r and q i the values of the statistic being 490 evaluated. For each AUC, we calculated a confidence interval by generating 10,000 bootstrap samples of 491 5,000 variants (with the same ratio of IBD:recurrent variants as the simulated sample). We then sorted the 492 10,000 AUC estimates and took the 2.5th and 97.5th percentiles to get a 95% confidence interval. 493 494

Calculating an expected fraction of recurrent mutations from polymorphism probabilities 495
As a proxy for the mutation rate, we estimated the polymorphism probability for 5-mer sequence contexts 496 (i.e., the focal base and two bases up and downstream) as the fraction of sites with that context that were 497 variable in the UK10K dataset. These polymorphism probabilities were highly correlated with those 498 calculated previously with the 1000 Genomes dataset (Aggarwala and Voight 2016) (Pearson's 499 correlation = 0.99, P<10 -100 ), with a higher fraction of polymorphic sites in the UK10K for each context 500 due to the larger sample size (Supplementary Figure 5). 501 To predict the fraction of sites that should be called recurrent based on sequence context 502 polymorphism probabilities, we used a simple Poisson model of mutation. With the polymorphism 503 probability for a context as the Poisson rate parameter λ and the number of mutations at a site H, the 504 probability of a recurrent mutation is the probability of two or more mutations at a site: 505 506 As we are only considering sites where there has been at least one mutation event, i.e. polymorphic sites, 507 the probability of a recurrent mutation at a site is then: 508

509
We calculated this probability for each 5-mer sequence context. We then calculated the expected fraction 510 by scaling the overall fraction of sites called non-IBD by each context's probability of a recurrent 511 mutation, relative to all the other contexts. 512 We used 5-mer sequence contexts for this analysis so that we would have a reasonable number of 513 variants classified as IBD or not for each sequence context. If we had used 7-mer sequence contexts, 514 some contexts would have too few variants to calculate the proportion called non-IBD. For the regression 515 models to predict non-IBD variants using multiple genomic annotations, we used 7-mer sequence 516 contexts, as there is significant mutation rate variation even within 5-mer contexts (Aggarwala and Voight 517 2016). 518 519

Identifying putative gene conversions 520
Within the set of variants called as non-IBD, we called putative gene conversion tracts that contained 2 or 521 more variants that were: 1) present in the same individuals, 2) at the same allele count, 3) within 1kb of 522 each other. 523 524

Genomic annotation datasets 525
We used the B statistic (McVicker et al. 2009) (downloaded from 526 http://www.phrap.org/othersoftware.html) to measure background selection, which estimates the 527 proportion of neutral variation in a region. VQSLOD and read depth were extracted from the UK10K 528 VCF files. We used a recombination rate map estimated for Europeans from the 1000 Genomes Project, 529 downloaded from 530 http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20130507_omni_recombination_rates (The 531 We can then use the binomial probability of the number of pairs with recombination distance greater than or less than the threshold distance to calculate the likelihood of the data given an IBD or recurrent allele.
2 BAYESIAN HIERARCHICAL MODEL TO IDENTIFY RECURRENT

MUTATIONS
As an alternative to the likelihood approach presented in the previous section, I also developed a Bayesian hierarchical model to jointly infer the model parameters along with identifying recurrent mutations. The benefit of this approach is that we can jointly infer the parameters of the posterior distribution and the posterior probability of each variant being recurrent from the data itself.
The model parameters are as follows: d : recombination distances for a variant t : TMRCA(s) for a variant (1 if IBD variant, 2-3 if recurrent depending on partition) k: component assignment of variant. k=1 if IBD, k>1 otherwise (each recurrent partition given its own value of k). j : given k, assignment of allele pairs. k = 1 if an IBD pair (k=1 for all pairs if k=1), k = 2 if a recurrent allele pair. µ, ø: mean and precision parameters for distribution of t . º: mixture fractions of each component k.
We can then use the binomial probability of the number of pairs with recombination distance greater than or less than the threshold distance to calculate the likelihood of the data given an IBD or recurrent allele.

MUTATIONS
As an alternative to the likelihood approach presented in the previous section, I also developed a Bayesian hierarchical model to jointly infer the model parameters along with identifying recurrent mutations. The benefit of this approach is that we can jointly infer the parameters of the posterior distribution and the posterior probability of each variant being recurrent from the data itself.
The model parameters are as follows: d : recombination distances for a variant t : TMRCA(s) for a variant (1 if IBD variant, 2-3 if recurrent depending on partition) k: component assignment of variant. k=1 if IBD, k>1 otherwise (each recurrent partition given its own value of k). j : given k, assignment of allele pairs. k = 1 if an IBD pair (k=1 for all pairs if k=1), k = 2 if a recurrent allele pair. µ, ø: mean and precision parameters for distribution of t . º: mixture fractions of each component k.

Rescaling the SFS with non-IBD mutations 538
Starting with the SFS calculated from all UK10K biallelic sites included in our study, for allele counts 2-5 539 for CpG->T and all other mutation types we calculated the fraction called as non-IBD. We then divided 540 this fraction by the power of our method, estimated by the percent of multiallelic sites identified as non-541 IBD at the chosen posterior threshold. From this fraction of non-IBD sites for the two mutation types, we 542 apportioned the non-IBD mutations into lower allele counts based on the relative frequency of allele 543 counts 1-5. For example, to determine what fraction of non-IBD 4-ton variants would be assigned 544 partition 1:3 vs. 2:2, we used the relative frequencies: 545 Where !:! is the relative frequency of the 1:3 partition, and ! is the frequency of singletons in the 547 original SFS. In the rescaled SFS, the number of singletons increased by the number of variants of allele 548 count 2-5 that were identified with partition 1:(n-1); the number of doubletons decreased by the number 549 of doubletons that were identified as recurrent, and increased by the number of variants of allele count 3-5 550 that had partition 2:(n-2); and so on through allele count 4. Allele count 5 was excluded from the rescaled 551 SFS plots because we did not identify recurrent variants at allele counts greater than 5.