Detection and quantification of inbreeding depression for complex traits from SNP data
Edited by Andrew G. Clark, Cornell University, Ithaca, NY, and approved June 21, 2017 (received for review December 22, 2016)
Letter
February 21, 2018
Letter
February 21, 2018
Significance
Inbreeding depression (ID) is the reduction of fitness in offspring of related parents. This phenomenon can be quantified from SNP data through a number of measures of inbreeding. Our study addresses two key questions. How accurate are the different methods to estimate ID? And how and why should investigators choose among the multiple inbreeding measures to detect and quantify ID? Here, we compare the behaviors of ID estimates from three commonly used SNP-based measures of inbreeding and provide both theoretical and empirical arguments to answer these questions. Our work illustrates how to analyze SNP data efficiently to detect and quantify ID, across species and traits.
Abstract
Quantifying the effects of inbreeding is critical to characterizing the genetic architecture of complex traits. This study highlights through theory and simulations the strengths and shortcomings of three SNP-based inbreeding measures commonly used to estimate inbreeding depression (ID). We demonstrate that heterogeneity in linkage disequilibrium (LD) between causal variants and SNPs biases ID estimates, and we develop an approach to correct this bias using LD and minor allele frequency stratified inference (LDMS). We quantified ID in 25 traits measured in participants of the UK Biobank, using LDMS, and confirmed previously published ID for 4 traits. We find unique evidence of ID for handgrip strength, waist/hip ratio, and visual and auditory acuity (ID between −2.3 and −5.2 phenotypic SDs for complete inbreeding; ). Our results illustrate that a careful choice of the measure of inbreeding combined with LDMS stratification improves both detection and quantification of ID using SNP data.
Sign up for PNAS alerts.
Get alerts for new articles, or get an alert when an article is cited.
Mating between close relatives has detrimental consequences on the survival and fertility of resulting offspring (1). This overall reduction of fitness, referred to as inbreeding depression (ID), is observable in a wide range of organisms, including plants (2), animals (3, 4), and humans (5). In humans, major abnormalities are more frequent in children from consanguineous marriages (6) and genes causing rare diseases can be mapped by ascertaining children from such matings (7). To date, although the genetic basis of ID is not completely elucidated, two main hypotheses are proposed to explain this phenomenon: homozygosity for partially recessive deleterious mutations and heterozygous advantage (overdominance) (1, 8). More generally, ID can be estimated for any complex trait, even if the trait is not an obvious component of fitness. For polygenic traits, ID can be detected if there is directional dominance (DD) across loci, which means that the phenotype of individuals who are heterozygous deviates from the average phenotypes of homozygous individuals in a consistent direction. For fitness components, DD is negative; i.e., on average homozygosity reduces fitness.
In practice, ID can be estimated from pedigree studies when the relationships between parents are known (6, 9). However, given the limited number and the small sizes of such studies in humans, contemporary efforts (5, 10) to quantify ID have instead used SNP genotyping platforms to directly estimate inbreeding coefficients (F). SNP data may allow a more accurate evaluation of inbreeding (11), in particular for distant and cryptic inbreeding, and allow inference to be drawn from large population data (10). Conceptually, once a measure of F is derived from SNP data, ID can subsequently be estimated by correlating phenotype with the estimated F.
Genome-wide estimators of F fall in two categories: average homozygosity measures across loci (irrespective of position) and measures of continuous runs of homozygosity (ROH). Using ROH, ID has been reported for diseases (12, 13), height (5), and cognition (10). ROH-based estimates of F () have been previously shown to better correlate with the unobserved pedigree-inbreeding coefficient compared with other measures of inbreeding (14, 15), which has made them a gold standard. However, the sampling variance of these estimates is large, and consequently large sample sizes (10) are required to detect ID with measures. In addition, estimation depends on arbitrary (although optimized) choices of multiple parameters like the minimum number of SNPs covered by a ROH, the distance between two consecutive ROHs, and the number of heterozygous genotypes allowed in each ROH. Setting ROH length cutoffs ignores the contribution to ID of smaller identity by descent segments due to distant ancestors.
Therefore, quantifying the theoretical properties (bias and variance) of ID estimates derived from is challenging. These two critical limitations led us to consider two other commonly used measures of inbreeding (3, 15), namely the excess of homozygosity inbreeding coefficient (hereafter denoted ) as estimated in PLINK (16) and the correlation between uniting gametes (hereafter denoted ) previously introduced as in Yang et al. (17), as potential efficient measures for detecting and quantifying ID. We present the theory underlying unbiased estimation of ID and compare through simulations the performances of these three measures of inbreeding. We then quantify ID in 25 quantitative traits measured in a large dataset of individuals from the UK Biobank, using an approach that is robust to different assumptions on the distribution of effect sizes, to possible directional effects of minor alleles and to population stratification.
Results
Theoretical Determinants of Unbiased Estimation of ID.
We assume that the phenotype of interest is a quantitative trait y with genetic component that is underlain by random additive and dominance effects of m independent causal variants. We denote as the expected depression in y resulting from complete inbreeding, where is the minor allele frequency (MAF) of the jth causal variant and the expectation of its dominance effect. In the absence of epistasis, fitness-related phenotypes linearly decrease with increasing inbreeding (Eq. S2). This well-established linear relationship naturally implies the use of linear regression methods to estimate ID.
Least-squares estimates of ID obtained with converge with increasing sample size toward . When explicitly calculating and with respect to the genotypes and the effect size distributions, we found under classical assumptions (Eq. S4 and Table S1) that is unbiased when the average linkage disequilibrium (LD) among observed SNPs equals the weighted (by effect sizes) average LD between causal variants and observed SNPs. Although influenced by the effect size distribution, the consistency of toward b is mainly driven by differences in LD between causal variants and observed SNPs. Therefore, a simple condition under which is unbiased is if the causal variants are a random subset of the observed SNPs. However, if the causal variants are enriched in high-LD regions of the genome, will overestimate the actual inbreeding depression. In contrast, if the causal variants are enriched in a low-LD region like DNAse-I hypersensitive sites or enhancers (18), or if they are enriched among low-frequency variants, is expected to underestimate the true effect. This is further illustrated in our first simulation (Fig. 1).
Fig. 1.
Table S1.
Alleles | aj | Aj | Allele frequency |
---|---|---|---|
ak | pjpk + Djk | (1 − pj) pk − Djk | pk |
Ak | (1 − pk) pj − Djk | (1 − pj) (1 − pk) + Djk | 1 − pk |
pj | 1 − pj | 1 |
aj and ak are the minor alleles of SNP j and SNP k, and Aj and Ak are their major alleles, respectively The MAFs of SNP j and SNP k are pj and pk, respectively. Djk is the LD parameter.
LD heterogeneity between causal variants and SNPs used for inference has been previously shown to determine the consistency of heritability estimates (19–21). We leveraged this estimation problem to propose a strategy to correct the differential LD bias when estimating ID. Following a previous approach (19), we explored how stratifying SNPs according to their LD score (22) and their MAFs before analyses (details given in Supporting Information) could correct or at least reduce these biases. We illustrate in our first simulation that LD score and MAF (LDMS) stratification performs well in correcting these biases (Fig. 1).
Similar to that shown with , we prove that the consistency of ID estimates obtained with (hereafter denoted ) is also determined by LD differences between SNPs and causal variants (Eq. S5). However, the bias of cannot simply be predicted by the ratio of the mean LD score in causal variants over the mean LD score in SNPs (Supporting Information). Nevertheless, our derivations predict that behaves similarly to with respect to LD differences between causal variants and SNPs. Importantly, we also prove that possible directional effects of minor alleles (DEMA) could confound because of the correlation between minor allele counts and (Supporting Information). Such directional effects could arise as a consequence of directional selection (when the minor allele is also the derived allele) as previously reported in human height (23) or simply because of population stratification (PS).
Simulation Study.
The complete description of the simulation study is given in Supporting Information.
Influence of differential MAF and LD between causal variants and SNPs.
We first considered three scenarios to illustrate the influence of LD and MAF heterogeneity between causal variants and SNPs. In all these scenarios, we assumed no DEMA, i.e., parameter in Eq. S3, and that b = −3 phenotypic SDs. Moreover, we assumed the expectation of the dominance effects to be either constant, i.e., , or inversely proportional to the variance of the minor allele count, i.e., . The first assumption corresponds to neutral traits, whereas the second one assigns a larger effect to SNPs with lower MAF and therefore corresponds more to traits under directional selection. Unbiasedness is defined below as when the average estimate of ID over multiple simulation replicates does not significantly differ from the value of b used for simulation.
Scenario 1.
In this scenario the causal variants were randomly sampled from the 3,857,369 autosomal SNPs that passed the genotypes quality control (Supporting Information). As predicted by our derivations, we observed that -based estimates of b were unbiased when dominance effects are assumed inversely proportional to the variances of allele counts, whereas an overestimation of of b was observed when dominance effects are assumed constant (Fig. 1A). This overestimation is explained by the fact that assuming a constant dominance effect, regardless of allele frequencies, creates an apparent MAF and LD heterogeneity between SNPs and causal variants by relatively up-weighting common SNPs compared with rarer SNPs (Eq. S3). We observed that LDMS stratification, which accounts for that heterogeneity, completely corrected this upward bias as presented in Fig. 1A. In addition, we found that produced unbiased estimates of b when dominance effects are assumed constant as for a neutral trait (Fig. 1A), but was biased downward ( of b) when dominance effects are inversely proportional to the variances of allele counts (Fig. 1B). This downward bias can be explained using the same reasoning presented above because in that case assuming dominance effects inversely proportional to the variances of allele counts relatively up-weights rarer SNPs compared with common ones. This downward bias could similarly be corrected using LDMS stratification. We also found that estimates of b obtained with ROH-based measures of inbreeding were strongly biased: of b using the definition of ROH from Joshi et al. (10) [hereafter denoted ] and of b using an alternative definition from Gazal et al. (15) or Howrigan et al. (24) [hereafter denoted as ]. The main difference between those two definitions of ROH is that requires LD pruning of the SNPs before calling the ROHs, whereas explicitly imposes a constraint on the ROH lengths (here Mb). This result highlights that LD pruning improves ID estimation using ROH-based inbreeding measures but still remains insufficient to produce unbiased estimates. Indeed, using more stringent LD pruning thresholds did not change our conclusion (Fig. S1). Overall, we found that LDMS stratified estimates for and were unbiased in all cases (Fig. 1 A and B), which emphasizes that this strategy can be safely used even when causal variants are perfectly tagged by SNPs.
Fig. S1.
On average over 1,000 simulation replicates we found that -associated estimates had smaller standard errors (SE) compared with or ( and ) (Fig. S2 A and B). consequently yielded the largest statistical power whereas was second best with a power on average below that of . On the other hand, because of their large SEs, and yielded the smallest statistical power to detect ID. Finally, we found that LDMS stratified estimates had larger SEs compared with nonstratified estimates. This increase in SE corresponds on average over all inbreeding measures to an loss of statistical power and is explained by the larger underlying effective dimensionality (4 LD score strata 6 MAF strata = 24 parameters actually estimated; Supporting Information) of the LDMS approach compared with the nonstratified inference.
Fig. S2.
Scenarios 2 and 3.
For the two other scenarios we used 1,358,699 SNPs within exons, introns, 3’-UTRs, 5’-UTRs, and promoter regions ±500 bp (SI Materials and Methods, URLs). SNPs within these five genomic (sets of) regions have distinct MAF and LD distributions as shown in Figs. S3 and S4. In scenario 2, we sampled the causal variants among 542,379 intronic SNPs with MAF whereas in scenario 3 causal variants were sampled among 28,341 SNPs within exons, 3’-UTRs, and 5’-UTRs. Our theoretical derivations predict that and would underestimate the true ID in scenario 2 because causal variants in that scenario had on average lower LD scores (Fig. S3). Accordingly, we found over 1,000 simulation replicates an underestimation of of the true ID for and (Fig. 1). These downward biases could be reduced below of b using LDMS stratification (Fig. 1 A and B) and were not significantly different from 0 . In scenario 3 causal variants had on average larger LD scores and MAF (Figs. S3 and S4). We therefore expected an overestimation of ID estimates in that scenario according to our theoretical derivations. This predicted upward bias was indeed more noticeable in our simulations ( on average over all inbreeding measures) compared with scenario 2. Still, using LDMS stratification we were able to reduce these biases down to of b on average over all inbreeding measures (Fig. 1) and more specifically of b for . Overall, LDMS stratification using yielded the smallest biases compared with all other strategies.
Fig. S3.
Fig. S4.
Influence of DEMA.
Let denote the expectation of the additive effect of the minor allele at the jth causal variant. We define as an overall measure of DEMA. Under this assumption we prove that estimates of ID obtained with are confounded because of the correlation between and minor allele counts (Supporting Information). We illustrate and quantify here that confounding bias. The parameters of this simulation are similar to scenario 1 of the first simulation, except that data are now simulated assuming no ID, i.e., and using . With , DEMA (which contributes to the additive genetic variance) accounts in our simulations for of the total phenotypic variance (Supporting Information). We considered two alternatives for the expectation of the additive effects : (i) is constant and (ii) is inversely proportional to the MAF . Under both alternatives, we found that estimates of ID obtained with were severely biased (between −2 and −3 phenotypic SDs whereas the true value is 0; Fig. 2A) unlike those derived from other inbreeding measures. Whether this bias can be corrected is a difficult question in practice. Indeed, under the simplistic scenario considered here where causal variants are well tagged, adjusting for a genetic score summing all minor alleles should completely correct this bias. However, in more realistic situations where causal variants are only partially tagged, this would remain insufficient. In contrast, theory shows that is orthogonal to minor allele counts and consequently would not be influenced by such directional effects if these exist. This result has motivated our decision to restrict real data analyses to only even if had better statistical power in our previous simulations.
Fig. 2.
DEMA could also arise as a consequence of PS. To illustrate that second aspect, we performed another simulation with settings similar to scenario 1 of the first simulation (with and ) but now include the effect of 10 genotypic principal components (PC) as a proxy for PS (Supporting information). When varying the contribution of PS to the phenotypic variance from 0 to , we observed that PS had a larger influence on estimates derived from and . These observations are consistent with our theoretical results (at least for ) as they directly derive from the correlation between and minor allele counts, the weighted sum of which constitutes the PC. The biases of - and -associated ID estimates were on average phenotypic SD (Fig. 2B) when PS explained of the phenotypic variance but were significant only for (Fig. S5B). On the contrary, the biases observed when using never exceeded 0.5 phenotypic SD (Fig. 2B) and were not statistically different from 0 (Fig. S5B). In conclusion, orthogonality between and minor allele counts (25) (Supporting Information) guarantees that confounding by DEMA or PS is negligible for ID estimates obtained with this measure.
Fig. S5.
To summarize this section, we show in our simulations that biases in ID estimates induced by MAF and LD heterogeneity between SNPs and causal variants can be corrected using LDMS stratification. Moreover, we show on average that LDMS correction performs better when applied to and that -based ID estimates are robust to DEMA and PS. Overall, offers the best trade-off between statistical efficiency and unbiasedness in the situations covered in this simulation study. We therefore recommend its use and focus hereafter our analyses on real data to .
Analysis of UK Biobank Data.
We quantified ID in 25 quantitative traits measured in (Table 1) participants from the UK Biobank (Supporting Information). These traits can be grouped into three categories including physical measures (standing height, weight, body mass index, waist and hip circumferences, waist/hip ratio, bone mineral density, body fat percentage, hand grip strength, systolic blood pressure, diastolic blood pressure, heart pulse rate, peak expiratory flow, visual acuity measured on log minimum angle of resolution (MAR) scale, and auditory acuity measured as the speech reception threshold), cognitive traits and educational attainment (fluid intelligence score, mean time to identify matches, maximum number of digits remembered, and age when completed full education), and sex-specific reproductive traits (number of children fathered, number of live births, age at menarche, and age at menopause) (Table S2). Some of these traits like standing height, peak expiratory flow (strongly correlated with forced expiratory volume in 1 s: , ), educational attainment, and cognitive ability were previously reported to be associated with inbreeding (10). Beyond quantifying the effect of inbreeding on these traits, we also aimed to evaluate whether differential LD and MAF distribution in causal variants influenced classical least-squares estimates and if so to correct these biases using our LDMS inference. The analyses were performed using linear regression adjusted for age, sex (for traits not specific to males or females), recruitment center, Townsend deprivation index (26) as a proxy for socio-economical status, and the first 10 genotypic PCs. The last three adjustments were considered to account for geographical and socio-economic structures in the UK population, which we found to correlate with levels of inbreeding (Table S3).
Table 1.
Traits | N | Estimate | SE | P value | Estimate | SE | P value | PHET |
---|---|---|---|---|---|---|---|---|
PEF | 117,575 (103,781) | −4.1 (−4.1) | 0.62 (0.74) | 4.57 (3.48 ) | −4.12 (−4.19) | 0.69 (0.84) | 2.25 (6.35 ) | 0.941 (0.81) |
AA, speech reception threshold | 43,175 (38,449) | 5.23 (4.44) | 1.04 (1.26) | 5.55 (4.45 ) | 5.34 (4.39) | 1.16 (1.45) | 4.6 (2.51 ) | 0.828 (0.94) |
FIS | 45,043 (40,089) | −3.9 (−3.14) | 0.97 (1.23) | 5.36 (1.08 ) | −4.72 (−4.32) | 1.06 (1.42) | 8.52 (2.25 ) | 0.061 (0.07) |
HGS, average of left and right hands | 139,623 (122,950) | −2.36 (−2.72) | 0.54 (0.68) | 1.38 (5.79 ) | −2.43 (−3.31) | 0.6 (0.76) | 4.64 (1.51 ) | 0.771 (0.1) |
NCF | 51,494 (45,483) | −3.09 (−3.69) | 0.96 (1.16) | 1.27 (1.44 ) | −4.01 (−4.58) | 1.07 (1.32) | 1.79 (5.17 ) | 0.039 (0.14) |
MTCIM | 138,902 (122,334) | 1.94 (1.7) | 0.54 (0.68) | 3.66 (1.21 ) | 2.05 (2.1) | 0.6 (0.77) | 6.14 (6.24 ) | 0.647 (0.26) |
WHR | 140,295 (123,540) | 2.35 (3.03) | 0.53 (0.67) | 1.12 (6.18 ) | 1.97 (2.89) | 0.6 (0.76) | 7.85 (1.37 ) | 0.121 (0.69) |
VA, log MAR scale | 29,616 (26,596) | 4.04 (5.77) | 1.20 (1.51) | 7.42 (1.37 ) | 4.39 (6.68) | 1.32 (1.73) | 9.04 (1.14 ) | 0.518 (0.26) |
Effect sizes and standard errors are expressed in phenotypic SD of the trait. Results presented in parentheses were obtained after removing 16,781 related individuals and 19 extreme cases of inbreeding (FUNI >0.15). N: sample size in the analysis. PHET is the P value from the LDMS heterogeneity test comparing nonstratified and LDMS-stratified estimates.
Table S2.
FUNI | FHOM | ||||||||
---|---|---|---|---|---|---|---|---|---|
Traits | N | Estimate (SE) | P value | Estimate (SE) | P value | Estimate (SE) | P value | Estimate (SE) | P value |
PEF | 117,575 | −4.1 (0.6) | 4.6 10−11 | −2.7 (0.5) | 2.1 10−8 | −7.6 (1.4) | 1.3 10−7 | −8.5 (1.5) | 3.9 10−8 |
AA | 43,175 | 5.2 (1.0) | 5.5 10−7 | 3.2 (0.8) | 7.5 10−5 | 9.4 (2.4) | 1.1 10−4 | 11 (2.6) | 1.3 10−5 |
WHR | 140,295 | 2.3 (0.5) | 1.1 10−5 | 1.8 (0.4) | 3.2 10−5 | 2.7 (1.2) | 0.02 | 4.5 (1.3) | 5.5 10−4 |
HGS | 139,623 | −2.4 (0.5) | 1.4 10−5 | −1.6 (0.4) | 1.5 10−4 | −3.8 (1.2) | 1.8 10−3 | −4.3 (1.3) | 1.1 10−3 |
FIS | 45,043 | −3.9 (1.0) | 5.4 10−5 | −2 (0.8) | 8.2 10−3 | −7.9 (2.2) | 3.9 10−4 | −7.9 (2.4) | 9 10−4 |
MTCIM | 138,902 | 1.9 (0.5) | 3.7 10−4 | 1.3 (0.4) | 2.4 10−3 | 3.7 (1.2) | 2.5 10−3 | 3.6 (1.3) | 6.2 10−3 |
Height | 140,302 | −1.9 (0.5) | 4.8 10−4 | −1.1 (0.4) | 0.01 | −2.7 (1.2) | 0.03 | −3.5 (1.3) | 7.3 10−3 |
VA | 29,616 | 4 (1.2) | 7.4 10−4 | 2 (0.9) | 0.04 | 7.7 (2.7) | 4.9 10−3 | 7.9 (2.9) | 7.3 10−3 |
NCF | 51,494 | −3.1 (1.0) | 1.3 10−3 | −2.4 (0.7) | 1.4 10−3 | −5.4 (2.2) | 0.01 | −7.7 (2.4) | 1.3 10−3 |
WC | 140,233 | 1.4 (0.5) | 7.2 10−3 | 0.91 (0.4) | 0.03 | 1.9 (1.2) | 0.12 | 2.7 (1.3) | 0.04 |
BFat | 138,055 | 1.1 (0.4) | 0.05 | 0.71 (0.4) | 0.09 | 1.6 (1.2) | 0.19 | 1.9 (1.3) | 0.14 |
HR | 130,694 | 1.1 (0.6) | 0.05 | 1.2 (0.4) | 5.7 10−3 | 2.6 (1.3) | 0.04 | 3.2 (1.4) | 0.02 |
MNDR | 145,75 | −3.4 (1.8) | 0.06 | −2.7 (1.4) | 0.06 | −6.7 (4.3) | 0.19 | −8.3 (4.5) | 0.07 |
A-Mena | 72,196 | 1.3 (0.8) | 0.08 | 0.29 (0.6) | 0.63 | 1.4 (1.7) | 0.42 | 1.9 (1.9) | 0.3 |
BMI | 139,686 | 0.94 (0.5) | 0.08 | 0.54 (0.4) | 0.21 | 0.86 (1.2) | 0.48 | 1.3 (1.3) | 0.3 |
BMR | 138,047 | −0.83 (0.5) | 0.12 | −0.61 (0.4) | 0.16 | −2.3 (1.2) | 0.05 | −2.5 (1.3) | 0.05 |
A-Edu | 95,540 | −0.6 (0.6) | 0.34 | −0.36 (0.5) | 0.48 | 0.19 (1.4) | 0.89 | −0.07 (1.5) | 0.96 |
NLB | 59,943 | 0.8 (0.9) | 0.36 | 0.26 (0.7) | 0.71 | 1.1 (2) | 0.57 | 0.77 (2.1) | 0.72 |
NIM | 138,543 | −0.45 (0.5) | 0.41 | −0.79 (0.4) | 0.07 | −2.9 (1.2) | 0.02 | −3.1 (1.3) | 0.02 |
A-Meno | 42,782 | −0.82 (1.0) | 0.43 | −0.2 (0.9) | 0.82 | −1.8 (2.5) | 0.46 | −2 (2.7) | 0.44 |
SBP | 130,787 | 0.32 (0.6) | 0.57 | 0.01 (0.4) | 0.99 | 0.24 (1.3) | 0.85 | 0.79 (1.4) | 0.56 |
BMD | 80,698 | −0.31 (0.7) | 0.65 | −0.26 (0.5) | 0.64 | 0.69 (1.5) | 0.65 | 0.08 (1.6) | 0.96 |
DBP | 130,831 | 0.17 (0.56) | 0.75 | −0.21 (0.4) | 0.63 | 0.8 (1.3) | 0.52 | 1.1 (1.4) | 0.39 |
HC | 139,694 | 0.11 (0.54) | 0.83 | −0.1 (0.4) | 0.81 | 0.5 (1.2) | 0.68 | 0.17 (1.3) | 0.89 |
Weight | 139,916 | −0.07 (0.5) | 0.89 | −0.096 (0.4) | 0.82 | −0.71 (1.2) | 0.56 | −0.68 (1.3) | 0.59 |
AA, auditory acuity measured as the speech reception threshold; A-Edu, age when completed full education; A-mena, age at menarche; A-Meno, age at menopause; BFat, body fat percentage; BMD, bone mineral density; BMI, body mass index; BMR, basal metabolic rate; DBP, diastolic blood pressure; FIS, fluid intelligence score; HC, hip circumference and weight; HGS, hand grip strength; HR, heart pulse rate; MNDR, maximum number of digits remembered; MTCIM, mean time to correctly identify matches during cognitive test; NCF, number of children fathered; NIM, number of incorrect matches during cognitive test; NLB, number of live births; PEF, peak expiratory flow; SBP, systolic blood pressure; VA, visual acuity in log MAR scale; WC, waist circumference; WHR, ratio of waist circumference over hip circumference. Effect sizes and SEs are expressed in phenotypic SD. N: sample size in the analysis. Other abbreviations are defined in Fig. S2. P values <0.05/25 = 0.002 are highlighted in boldface type.
Table S3.
Measure of population structure | F statistic | Degrees of freedom | P value |
---|---|---|---|
Age at recruitment | 6.84 | 1, N | 8.9 10−3 |
Sex | 0.32 | 1, N | 0.57 |
Recruitment center | 2.68 | 21, N | 4.5 10−5 |
Townsend deprivation index Genotypic PCs | 47.6 | 1, N | 5.4 10−12 |
PC1 | 94.6 | 1, N | 2.4 10−22 |
PC2 | 202 | 1, N | 7.9 10−46 |
PC3 | 19.4 | 1, N | 1.1 10−5 |
PC4 | 6.53 | 1, N | 0.01 |
PC5 | 1.66 | 1, N | 0.19 |
PC6 | 3.15 | 1, N | 0.07 |
PC7 | 0.71 | 1, N | 0.40 |
PC8 | >3.35 | 1, N | 0.07 |
PC9 | 0.03 | 1, N | 0.86 |
PC10 | 1.63 | 1, N | 0.20 |
Statistical association was measured analysis of variances. N = 140,720.
After Bonferroni correction ( traits = ), we detected significant ID in eight traits (Table 1), using LDMS stratified inference based on . Ranked by decreasing magnitude of depression, these traits are auditory acuity (AA), fluid intelligence score (FIS), visual acuity (VA), peak expiratory flow (PEF), number of children fathered (NCF), hand grip strength (HGS), mean time to correctly identify matches (MTCIM), and waist/hip ratio (WHR). This analysis included 16,781 related individuals (estimated to be first, second, or third degree) and 19 participants with extreme inbreeding . As a sensitivity analysis we reran all analyses without related and inbreeding outliers, which reduced the number of traits passing the significance threshold to five. The three traits dropped in this secondary analysis were AA ( in secondary analysis), FIS ( in secondary analysis), and MTCIM ( in secondary analysis). To test whether the differences in ID estimates between full and reduced analyses are significant, we used a jackknife procedure to compare the observed differences with differences generated when randomly excluding 16,800 participants. Over 1,000 resampling events, we found that the observed differences in ID estimates for the eight traits highlighted above were not significantly different from those obtained when excluding random subsets (empirical ; Fig. S6). We consequently believe that the drop of significance between those two analyses is mainly explained by the reduced statistical power and not by confounding. In addition, we explored how much of ID could be captured at genome-wide significant (GWS) SNPs. We therefore selected trait-specific GWS SNPs from the genome-wide association studies (GWAS) catalog (SI Materials and Methods, URLs) and assessed inbreeding depression for the same traits using at those loci. We could not, however, detect any significant association with the traits analyzed in our study. Even for height, for which common GWSs are now reported (27), the estimate of inbreeding depression at GWS was only −0.08 SD for complete inbreeding .
Fig. S6.
We observed for all traits that ID estimates derived from were systematically larger than those obtained with (Table S3). As expected, their SEs were also larger. In particular, only four and six traits (of eight detected with ) passed the Bonferroni threshold when using and , respectively. On the other hand, ID estimates obtained with were systematically smaller than those obtained using (Table S3), with an average over the eight traits significantly associated with . The latter observation would be expected if the traits analyzed are under directional selection as observed in our simulations when rarer variants were assumed to have larger effects.
We observed for most traits that LDMS stratified and nonstratified estimates were similar (Table 1), suggesting weak differential LD and MAF distributions in SNPs tagging causal variants. Nonetheless, a marginally significant (LDMS heterogeneity test ; Table 1 and Fig. S7) difference could be observed in NCF for which the LDMS ID estimate was SE larger than the nonstratified one (Table 1). This also translated into an improvement of the association P value from to (Table 1). We subsequently assessed which component(s) in the LDMS stratification contributed the most to NCF (Fig. S8). We therefore fitted a first multivariate regression model adjusted for four inbreeding coefficients specific to each LD score strata component and then another multivariate regression model adjusted for six inbreeding coefficients specific to each MAF stratum. We chose to fit two different models (for MAF and LD separately) instead of one including 24 covariates to minimize the effects of colinearity between inbreeding measures in each LDMS stratum. We found a nominally significant contribution of SNPs with minor alleles ( = −4.01 phenotypic SD; ) but no significant enrichment in LD strata despite a large contribution of the second-lowest LD score stratum ( = −1.62 phenotypic SD; ). According to our derivations, this enrichment of ID in lower-frequency SNPs and more generally in low-LD regions explains why nonstratified analyses produced smaller estimates compared with the LDMS approach. These results imply a disproportionate contribution of low-frequency SNPs to ID in NCF.
Fig. S7.
Fig. S8.
SI Materials and Methods
Statistical Models and Notations.
We consider the following model:For individual , is the observed value of the fitness component of interest, and is the minor allele count at the th causal SNP (). is the MAF of the th causal SNP, is the indicator of heterozygosity, and is a residual term capturing nongenetic effects on the observed phenotype. The additive and dominance effect sizes of the minor allele at the th causal SNP are denoted and , respectively. For the sake of simplicity, we assume that SNP interactions are negligible and consequently do not appear in Eq. S1. We also assume independence between the causal variants, between the genotypes and the effect sizes, and between the genetic and the nongenetic effects. Finally we assume the effect sizes to be random and such that and . We denote individual s inbreeding coefficient. From Eq. S1 the expectation of with respect to the genotypes and the effect sizes can be written aswhich translates a linear and often negative relationship between the average fitness component and the inbreeding coefficient . The slope in this linear relationship is the expected depression in the fitness component due to complete inbreeding (). Eq. S1 can also be written aswhere and are the nondirectional additive and dominance breeding values of individual , respectively. and are defined asand and are defined as
[S1]
[S2]
[S3]
Measures of Inbreeding.
We studied three measures of inbreeding. All these measures of inbreeding require individual SNP genotypes and can be used in the absence of any pedigree information. The first inbreeding measure is the excess of homozygosity measure defined here aswhere is the minor allele count of SNP , is the MAF, and is the number of genotyped or imputed SNPs available. is implemented in PLINK2 (command: –het).
The second measure () is based on the correlation between uniting gametes. This measure was defined in Yang et al. (17) as is implemented in PLINK2 and GCTA software (command: –ibc). The last measure is defined as the proportion of the genome within ROH. More specifically, the inbreeding measure was calculated as the cumulated length (in base pairs) of one individual’s genome within ROHs divided by (the approximate length of the autosomal genome in base pairs). We used two definitions of ROHs corresponding to those proposed in Joshi et al. (10) (definition 1) and Gazal et al. (15) (definition 2). Inbreeding measures calculated using definition 1 and definition 2 are respectively denoted as and . Both definitions used SNPs with MAF % as previously reported in Joshi et al. (10) and Gazal et al. (15).
Definition 1 requires ROHs to be at least 1.5 Mb long, to cover at least 50 SNPs, and to be located at least 1 Mb away from one another. This corresponds to the following PLINK2 command: –homozyg –homozyg-density 50 –homozyg-gap 1000 –homozyg-kb 1500 –homozyg-snp 50 –homozyg-window-het 1 –homozyg-window-missing 5 –homo- zyg-window-snp 50. Definition 2 requires only LD pruning the SNPs (PLINK2 command –indep 50 5 2), but does not impose additional explicit constraint on ROHs besides not allowing any heterozygote in the ROHs. The following PLINK1 command was used to call ROHs under definition 2: –homozyg –homozyg-window-het 0 –homozyg-snp 50 –homozyg-kb 0 –homozyg-density 5000 –homozyg-gap 5000. As explained in Gazal et al. (15), the options –homozyg-density 5000 and –homozyg-gap 5000 are set for PLINK1 to ignore those constraints.
Theoretical Determinants of Unbiased Estimation.
The linear relationship in Eq. S2 naturally implies the use of linear regression to estimate ID. Least-squares estimates of ID derived from the inbreeding measure F asymptotically converge toward .
For F = , can be written as (proof in SI Appendix)In Eq. S4, is the average LD between the th causal variants and the observed SNPs, whereas is the average LD between the th observed SNPs and the observed SNPs. To better understand Eq. S4, we note that when is inversely proportional to , simplifies into and therefore is simply the ratio of mean LD score in causal variants over the mean LD score in SNPs. More generally, if (which would hold when causal variants are a subset of SNPs), then Eq. S4 simplifies as .
[S4]
For F = we can derive a somewhat similar expression (proof in SI Appendix),In Eq. S5, and measure the LD between causal variant and SNPs, whereas measures the LD among SNPs. Detailed derivations of these quantities are given in SI Appendix as well as the proof of unbiasedness when causal variants are randomly sampled from SNPs.
[S5]
LD Score and MAF (LDMS) Stratified Inference.
We defined the LD score (22) of SNP j as the sum of pairwise between SNP and all SNPs located within 10 Mb around SNP ’s physical position. To correct potential biases induced by uneven LD (and MAF) distribution between causal variants and SNPs, we used a LD and MAF stratified inference as previously described (19). This approach is summarized below in five steps.
(i) The first step consists of identifying segments of the genome with uniform LD patterns. For this we used a sliding-window approach where each segment contains SNPs with overlap of SNPs with adjacent segments. The number was chosen as in Yang et al. (19) to be the average number of SNPs per 100 kb of a chromosome: , with and being respectively the number of SNPs and the length (in kilobases) of the chromosome. We then calculated the mean LD score of the variants in each segment and took the average of the mean LD scores of two adjacent segments for overlapped regions. This process partitioned the genome into a large number of segments. (ii) The second step consists of stratifying the SNPs into K groups according to the mean LD score of the segments () and the MAF interval they belong to. We used here four classes of LD [first quartile (low LD), first quartilemedian, median 3rd quartile, and third quartile) and six MAF intervals [(1–5%), (5–10%), (10–20%), (20–30%), (30–40%), and (40–50%)], i.e., strata in total. (iii) For each group , the third step consists of calculating a measure of inbreeding (could be any of the three measures presented in this study) specific to that group. (iv) In the fourth step we fitted a multivariate linear regression model using all of the s as covariates to obtain , a vector of regression coefficients corresponding to each group. (v) The fifth and last step comprises collapsing these estimates into a single genome-wide estimate : , where is a -dimensional vector with all entries equal to 1. The SE of ) can be obtained from the asymptotic variance–covariance matrix of , using the following relationship:We developed and applied a Wald test to assess the statistical significance of differences between nonstratified and LDMS stratified ID estimates. The test statistic is defined as the squared difference between the two estimates divided by the variance of that difference. The details of that derivation and some additional simulations showing that this test is not inflated in the absence of LD or MAF heterogeneity are presented in SI Appendix and Fig. S7. The LDMS procedure is implemented in GCTA as part of the GCTA-GREML-LDMS analysis (SI Materials and Methods, URLs).
UK Biobank Data.
Inclusion criteria.
We used baseline data from 152,729 participants from the UK Biobank (34). To ensure ancestry homogeneity, we selected individuals who reported to be British, Irish, white, or of any other white background and whose coordinates on the first genetic PC were below 0 (Fig. S9). In total, we included 140,720 participants in this analysis. All participants in the UK Biobank study provided written informed consent.
Phenotyping.
We analyzed 25 traits: standing height, weight, body mass index, waist and hip circumferences, ratio between waist circumference and hip circumference (referred to as WHR), bone mineral density, body fat percentage, HGS, systolic blood pressure, diastolic blood pressure, heart pulse rate, PEF, VA (measured on log MAR scale) and AA (measured as the speech reception threshold), cognitive traits and EA (FIS, MTCIM, maximum number of digits remembered, and age when completed full education), and sex-specific reproductive traits (NCF, number of live births, age at menarche, and age at menopause). For each trait, we excluded values SD away from the mean. The sample sizes specific to each trait are given in Table S2. Finally, all traits were adjusted for age, sex (for traits not specific to males or females), recruitment center, 10 genotypic PCs, and Townsend deprivation index (26) as a proxy for the socio-economic status. After adjustment, we used inverse normal transformed residuals as the traits to be correlated with inbreeding measures. Although focused on , we also provide ID estimates from other inbreeding measures, in particular from for comparison with previous and future studies using those measures (Table S2).
Genotyping.
We used data from 152,729 participants (including 140,720 analyzed in this study; SI Materials and Methods, Inclusion criteria) genotyped in the first phase of genotyping of the UK Biobank. The first steps of the quality control have been previously described (SI Materials and Methods, URLs). Phasing and imputation were performed using SHAPEIT and IMPUTE2 (SI Materials and Methods, URLs), respectively, as previously described (35). After imputation, we selected 9,493,148 autosomal SNPs with imputation quality , MAF1%, and Hardy–Weinberg equilibrium test P value . Imputed SNPs were then called to the genotypes having the largest posterior probability. Finally, we removed redundancy by LD pruning SNPs with a squared genotype correlation . In total we used 3,857,369 SNPs in this analysis.
Simulation Study.
In all simulations we used 10,000 unrelated (pairwise genetic relationship ) individuals randomly sampled from the 140,720 individuals selected for analysis.
Simulation 1: Influence of MAF and LD distribution in causal variants.
To illustrate the importance of LD and MAF distribution in causal variants, we considered three scenarios. In the first scenario (scenario 1) we sampled the causal variants at random among the 3,857,369 observed SNPs. For the two other scenarios we used 1,358,699 SNPs within exons, introns, 3-UTR, 5-UTR, and promoter regions 500 bp (SI Materials and Methods, URLs). The LD score and MAF distributions of SNPs within these functional annotations are presented in Figs. S3 and S4. In scenario 2, we sampled the causal variants among 542,379 intronic SNPs with MAF % (low-LD scenario) whereas in scenario 3 causal variants were sampled among 28,341 SNPs within exons and 3- and 5-UTRs (high-LD scenario). In all simulations, we sampled 1,000 sets of causal variants each and simulated a phenotype (1 replicate per set of causal variants), using the following modified version of Eq. S3:In Eq. S6, and are respectively defined as in ref. 25 as the scaled additive and orthogonal dominance values:Under Eq. S1 and assuming that and , the phenotypic variance can be expressed aswhere is the variance of the nongenetic effects [derivation of Eq. S7 in SI Appendix]. Zhu et al. (25) found on average across multiple traits that additive variance was about fivefold larger than dominance variance. Following their findings, we fixed in our simulations , , and , where denotes the heritability of the trait. We observed in our simulations that was negligible and consequently . We used two values of , and SD of the trait, to study statistical power to detect ID. In all simulations, the statistical significance threshold for P values is 5%. The ID parameter was estimated for each measure of inbreeding by fitting a linear regression model with the phenotype as the response variable and the corresponding measure as the primary covariate.
[S6]
[S7]
Simulation 2: Influence of directional effects of minor alleles.
In this simulation we illustrated how directional effects of minor alleles may confound estimates of ID. The settings of this simulation are similar to scenario 1 in the first simulation. We simulated data assuming no inbreeding depression, i.e., and fixed . We considered two alternatives for the expectation of the additive effects : (i) is constant and (ii) is inversely proportional to the MAF .
We further quantified the effect of PS on ID estimates by running another simulation (also similar to scenario 1 in the first simulation) with , , and now including the effect of the first 10 genotypic PCs as a proxy for PS:In Eq. S8 was chosen so that PS accounts for 1%, 2%, 3%, 4%, and 5% of the total phenotypic variance. PCs were calculated using genotypes from the 10,000 unrelated participants selected for the simulations. Type I error in this simulation was calculated as the fraction of P values (testing for the significance of ID) below 0.05.
[S8]
Selection of Genome-Wide Significant SNPs from the GWAS Catalog.
To assess how much ID could be captured at GWS SNPs, we assessed inbreeding using at GWS SNPs from the GWAS catalog (SI Materials and Methods, URLs). We downloaded version 1.0 of the GWAS catalog and selected SNPs so that their associated “DISEASE/TRAIT” field exactly matched “height,” “weight,” “body mass index,” “waist circumference,” “waist-hip ratio,” “bone mineral density,” “body fat percentage,” “systolic blood pressure,” “diastolic blood pressure,” “intelligence” [or “intelligence (childhood)”], “cognitive test performance,” “cognitive function,” “cognitive performance,” “general cognitive ability,” “educational attainment,” “educational attainment (years of education),” “educational attainment (college completion),” “age at first birth,” “menarche (age at onset),” “menopause (age at onset),” “forced expiratory volume in 1 s,” and “lung function (forced expiratory flow between 25% and 75% of forced vital capacity).” In total we selected 2,069 unique SNPs.
URLs.
i)
UK Biobank quality control: www.ukbiobank.ac.uk/wp-content/uploads/2014/04/UKBiobank_genotyping_QC_documentation-web.pdf.
ii)
UK Biobank imputation: www.ukbiobank.ac.uk/wp-content/uploads/2014/04/imputation_documentation_May2015.pdf.
iii)
SHAPEIT3, IMPUTE2: https://jmarchini.wordpress.com/software/.
iv)
Genomic annotations from UCSC website: https://genome.ucsc.edu/cgi-bin/hgTables.
vi)
GWAS catalog: https://www.ebi.ac.uk/gwas/.
SI Appendix
Derivation of Equation S2.
Derivation of Eqs. S4 and S5.
Additional notations and useful formulas.
We introduce now new notations. Let and be the number of minor alleles (0 or 1) inherited from the father and the mother, respectively. Then the minor allele count is defined as . We introduce as the LD parameter between minor alleles of SNP and SNP . We assume that the LD between paternal and maternal alleles is the same, i.e., that and have the same probability distribution. The latter distribution is presented in Table S1.
Let be the indicator of complete inbreeding; i.e., with probability if individual is completely inbred and if its parents are completely independent. We note . We have now thatMoreover, because we assume that LD between the father’s alleles is the same as between the mother’s alleles, then we have that .
Another important formula isFollowing the same reasoning, we can prove thatwhere is defined as
[S9]
Getting back to Eqs. S4 and S5.
Eq. S4.
Eq. S5.
For we have thatwhereTherefore, with , we have thatwhereand and both measure LD between the th causal variant and the SNPs used for inference, and measures LD between SNP and all observed SNPs. To better understand the relationship described in Eq. S12, let us consider a special case. When and is constant, Eq. S12 can be writtenwhere . Therefore, when causal variants are randomly sampled from observed SNPs, and , which guarantees unbiasedness.
[S11]
[S12]
[S13]
Derivation of Eq. S7.
To derive Eq. S7, we assumed the following modified version of Eq. S3:In Eq. S14 and are respectively defined as the scaled additive and orthogonal dominance values (25):Because of the independence between effect sizes and genotypes, we have thatSimilarly, we have that . Therefore,Let us now derive the expression of . We note that :We consider two cases corresponding to the different configurations presented in this paper:
[S14]
.When , we have thatIf we assume for example the causal variants to be a random subset of HapMap 3 (HM3) SNPs (1.1 million SNPs), we can therefore use their allele frequencies to approximate . Using allele frequencies of HM3 SNPs in the UK Biobank, we found that and therefore thatA different set of SNPs could have been used instead of HM3 SNPs. However, we argue here that in any case, is constant with respect to the number of causal variants and therefore would be negligible for a trait with a large number of causal variants. When , we have thatUsing HM3 SNPs, we can similarly approximate and , which givesIn conclusion, we show here that is of order and therefore contributes little to the phenotypic variance of traits with large numbers of causal variants.
.When , we have (still using HM3 SNPs) thatand similarly, when , we can write thatJust like , is also of order and hence contributes little to the phenotypic variance.
We now denote and asLet and be defined aswith . Hence,Now assuming independence between the causal variants we can derive Eq. S7,where is constant with respect to the number of causal variants. For large , .
Effect of Strong PS.
The contribution of PS reported in this study needs, however, to be put in perspective as our simulations and real data analyses were based upon a relatively homogeneous population within the United Kingdom. For stronger PS (e.g., between European and Asian populations) -based ID estimates can also be biased. This extreme situation would in general be handled as part of quality control. Still, we quantify below in a simulation study the bias of -based estimates in such a context.
Description of the simulation study.
Model and notations.
We assume that the data are collected from two populations: population 1 and population 2. For the sake of simplicity, we assume that all causal variants are observed and have the same allele frequency in population 1 and in population 2. We denote as and as the proportion of samples from population 2 in the data used for inference.
For individual in population ( or ), we simulated a trait from the modelwhere and is the MAF of SNP in population ; i.e., and . We defined as so that the additive genetic variance equals without directional effect because the sum of all is 0. We also defined , where is the depression for complete inbreeding. The residual term was simulated from a normal distribution with mean 0 and variance equal to , where is the heritability of the trait.
Simulation parameters.
In this simulation we fixed , , , and the total sample size . We varied between and and the proportion of samples from population 2 between 0.1 and 0.5.
Simulation results.
Table S4 reports the averaged estimates (over 10 simulation replicates) of inbreeding depression for the different combinations of and . We observed for large values of that produced biased estimates of inbreeding depression, but that bias was at most of the true effect. On the contrary, shows larger biases (on average 50% of ), which is in line with our results presented in the main text.
LDMS Heterogeneity Test.
Derivation of the test statistic.
We developed and applied a Wald test to assess the statistical significance of a difference between nonstratified ID estimation and LDMS stratified ID estimation.
Let us denote the matrix of inbreeding measure for each individual (–) and each SNP (–). The nonstratified measure of inbreeding (subscript means genome) is defined as the mean inbreeding across all SNPs. Therefore, , where is the -dimensional vector with all entries equal to 1.
Stratifying the SNPs into LD and MAF categories defines a linear map of the individual SNP inbreeding measures into aggregated inbreeding measures for each of the strata. The matrix of these strata-specific inbreeding measures can be writtenwhere is a matrix with entries so that equals the proportion of SNPs in the th category if SNP belongs to that category or equals 0 otherwise.
Without loss of generality, we assume that the sample means of and the columns of are 0. Using least-squares estimation theory, we can write the closed-form expressions of the the nonstratified () and the LDMS stratified () estimators. Assuming is the vector of the phenotypes (with sample mean 0), we have thatWe note that both and are linear functions of . We can therefore use the following matrix relationship:When the individuals in the study are independent and when assuming homoscedasticity of the residuals, we have that . We therefore deduce that the correlation between and isTherefore, to test the hypothesis vs. , we can use the test statisticwhere and are respectively the standard errors of and .
[S15]
[S16]
[S17]
Control of type I error and statistical power.
To check that type I error is well controlled under the null hypothesis, we analyzed the distribution of P values from the comparison of with when the causal variants are a random sample of the SNPs. Under this scenario, we expect and to coincide as illustrated in our first simulation study (simulation 1, scenario 1). We reused data generated in our first simulation study but simply varied the value of the true ID parameters between 0, –3, –6, and –9 SDs of the simulated trait. The results presented in Fig. S7 correspond to the comparison of with under the assumption that the dominance effect at each SNP is inversely proportional to the variance of allele counts.
As expected, the P-values distribution in this scenario was indistinguishable from the uniform distribution (Kolmogov–Smirnov test for all values of ). In particular, the observed proportion of P values below 5% does not deviate from its expectation of 5% (Fig. S7B). This result confirms that this test can be used safely to detect LD and MAF heterogeneity.
When we simulated LD and MAF heterogeneity by sampling causal variants from high or LD regions of the genome (simulation 1, scenarios 2 and 3), we noticed a deviation from the uniform distribution as illustrated in Fig. S7C. The power calculations presented in Fig. S7D are based on the fixed sample size of used in our simulations.
Discussion
We comprehensively quantified the behavior of ID estimators based on three commonly used measures of inbreeding. Our study illustrated some of the shortcomings of the most commonly used ROH-based estimates of ID, which not only are biased but also have large SE (approximately three times larger compared with ). This, along with the arbitrary choices underlying the definition of ROHs, leads us to recommend the use of over . Overall, our results suggest that -based ID estimates are robust to different assumptions about the distribution of effect sizes, to possible directional effects of minor alleles, and also to population stratification. The contribution of population stratification reported in this study needs, however, to be put in perspective as our simulations and real data analyses were based upon a relatively homogeneous population within the United Kingdom. For stronger population stratification (e.g., between European and Asian populations) -based ID estimates can also be biased. This somewhat extreme situation, which would in general be handled as part of quality control, is discussed in Supporting Information (Table S4 and Dataset S1).
Table S4.
= 10% | = 25% | = 50% | ||||
---|---|---|---|---|---|---|
Fst | FUNI (SE) | FHOM (SE) | FUNI (SE) | FHOM (SE) | FUNI (SE) | FHOM (SE) |
−0.10 | −0.88 (0.01) | −0.23 (0.005) | −0.87 (0.009) | −0.13 (0.003) | −0.93 (0.01) | −0.09 (0.004) |
−0.05 | −0.99 (0.01) | −0.59 (0.005) | −0.99 (0.006) | −0.43 (0.005) | −0.99 (0.01) | −0.35 (0.005) |
0.0 | −1.00 (0.01) | −1.00 (0.006) | −1.00 (0.008) | −1.00 (0.007) | −1.00 (0.01) | −1.00 (0.004) |
+0.05 | −0.98 (0.01) | −0.71 (0.007) | −0.99 (0.008) | −0.56 (0.006) | −0.99 (0.01) | −0.47 (0.004) |
+0.10 | −0.89 (0.01) | −0.45 (0.005) | −0.89 (0.01) | −0.30 (0.004) | −0.95 (0.01) | −0.24 (0.003) |
The true inbreeding depression parameter used to simulate the data is b ⁼ −1. Datasets are simulated as mixtures of two populations with mixture proportions π and 1 − π and allele frequency difference equal to Fst.
This study also highlights that differential LD distribution between causal variants and SNPs could bias ID estimates. As previously reported for genomic-relatedness-based restricted maximum-likelihood (GREML) (28) heritability estimates (19–21), we demonstrate through simulations that an LDMS approach successfully corrects these biases when estimating ID. More generally, the flexibility of the LDMS approach in terms of numbers and types of MAF/LD strata allows adaptation to any effect size distribution. Indeed, all SNP-based inbreeding measures are defined upon an underlying assumption on the distribution of dominance effects (e.g., when assuming constant dominance effects, the underlying inbreeding measure is ), which, when not verified, creates biases in ID estimates even when causal variants are randomly distributed among observed SNPs. We showed in our simulations using two distributions of dominance effects that such biases are explained by MAF and LD heterogeneity between causal variants and SNPs and therefore can be corrected using the LDMS approach. This result is important as it guarantees an unbiased estimation of ID regardless of the distribution of dominance effects.
Beyond methodological considerations, we confirmed in this study known associations between increased inbreeding and reduced lung function (PEF), cognitive ability (FIS and MTCIM), and fertility (NCF), which were previously reported, however, using different proxy traits (10, 29–31). We also replicated the association between inbreeding and decreased height (: −1.71 phenotypic SD for complete inbreeding; ) even though it was below the Bonferroni threshold. We did not replicate the association with educational attainment (EA) measured, as the “age when completed full education.” However, when measuring EA as whether or not participants went to college or university as in Joshi et al. (10), we found a strong although nominally significant negative association (odds ratio of 0.04, ).
We report evidence of ID for HGS, VA, AA and WHR. Although HGS, VA, and AA seem obvious fitness-related traits, these results still require replication. The association with WHR is particularly interesting as it illustrates on real data that one may benefit from using a less variable measure of inbreeding. Joshi et al. (10) reported a positive effect of inbreeding on WHR using ; however, even with individuals the effect did not reach statistical significance . Although heterogeneity between the cohorts involved in that meta-analysis may explain that apparent lack of statistical power, our theoretical and simulation results predict that if the authors had used a different metric [ instead of to reduce the SE] combined with a LD and MAF stratified inference, then such an effect would have been easily detected. We have not considered in this study the detection and estimation of nonlinear effects of inbreeding because of a lack of statistical power to detect such effects. Nonlinearity is predicted by theory in the presence of epistasis involving dominance (32) and was implied by Szpiech et al. (33), who showed that long runs of homozygosity were enriched for coding variants that are predicted to be deleterious. In conclusion, we have demonstrated that LD and MAF stratified inference based on as a measure of inbreeding minimizes bias relative to the other ID estimation strategies compared in this study. As illustrated here on real data, we believe that our approach will lead to more discoveries in forthcoming and larger studies.
Materials and Methods
Statistical Models and Notations.
We consider the following model:For individual i, is the observed value of the phenotype of interest, and is the minor allele count at the jth causal SNP (). We denote the minor allele frequency of the jth causal SNP, is the indicator of heterozygosity, and is a residual term capturing nongenetic effects on the observed phenotype. The additive and dominance effect sizes of the minor allele at the jth causal SNP are respectively denoted and . We assume independence between the m causal variants, between the genotypes and the effect sizes, and between the genetic and the nongenetic effects. Finally, we assume the effect sizes to be random and such that and .
[1]
Measures of Inbreeding.
We studied three measures of inbreeding. All these measures of inbreeding require individual SNP genotypes and can be used in the absence of any pedigree information. The first inbreeding measure is the excess of homozygosity measure defined here aswhere is the minor allele count of SNP k, is the minor allele frequency, and p is the number of genotyped or imputed SNPs available. is implemented in PLINK2 (command: –het).
The second measure () is based on the correlation between uniting gametes. This measure was defined in Yang et al. (17) as is implemented in PLINK2 and GCTA software (command: –ibc).
The last measure is defined as the proportion of the genome within ROH. More specifically, the inbreeding measure was calculated as the cumulated length (in base pairs) of one individual’s genome within ROHs divided by (the approximate length of the autosomal genome in base pairs). We used two definitions of ROHs corresponding to those proposed in Joshi et al. (10) (definition 1) and Gazal et al. (15) (definition 2). Inbreeding measures calculated using definition 1 and definition 2 are respectively denoted as and . Both definitions used SNPs with as previously reported in Joshi et al. (10) and Gazal et al. (15) (Supporting Information).
UK Biobank Data.
We used baseline data from 152,729 men and women who were genotyped in the first phase of genotyping of the UK Biobank (34). To ensure ancestry homogeneity, we selected individuals who reported to be “British,” “Irish,” “white,” or of “any other white background” and whose coordinates on the first genetic PC were below 0 (Fig. S9). In total, we included 140,720 participants in this analysis. The Northwest Multicentre Research Ethics Committee (MREC) approved the study and all participants in the UK Biobank study provided written informed consent. The first steps of the quality control have been previously described (SI Materials and Methods, URLs). Phasing and imputation were performed using SHAPEIT and IMPUTE2 (SI Materials and Methods, URLs), respectively, as previously described (35). After imputation, we selected 9,493,148 autosomal SNPs with imputation quality , , and Hardy–Weinberg equilibrium test P value . Imputed SNPs were then called to the genotypes having the largest posterior probability. Finally, we removed redundancy by LD pruning SNPs with a squared genotype correlation . In total we used 3,857,369 SNPs in this analysis.
Fig. S9.
Acknowledgments
This research was supported by the Australian Research Council (DP130102666 and DP160103860), the Australian National Health and Medical Research Council (1078037 and 1107258), the National Institute of Health (NIH Grants R01AG042568 and P01GM099568), and the Sylvia and Charles Viertel Charitable Foundation. This research has been conducted using the UK Biobank Resource under Project 12505.
Supporting Information
Supporting Information (PDF)
- Download
- 6.15 MB
Dataset_S01 (PDF)
- Download
- 62.27 KB
References
1
D Charlesworth, JH Willis, The genetics of inbreeding depression. Nat Rev Genet 10, 783–796 (2009).
2
X Huang, et al., Genomic analysis of hybrid rice varieties reveals numerous superior alleles that contribute to heterosis. Nat Commun 6, 6258 (2015).
3
J Huisman, LEB Kruuk, PA Ellis, T Clutton-Brock, JM Pemberton, Inbreeding depression across the lifespan in a wild mammal population. Proc Natl Acad Sci USA 113, 3585–3590 (2016).
4
JM Pemberton, PE Ellis, JG Pilkington, C Bérénos, Inbreeding depression by environment interactions in a free-living mammal population. Heredity 118, 64–77 (2017).
5
R McQuillan, et al., Evidence of inbreeding depression on human height. PLoS Genet 8, e1002655 (2012).
6
AH Bittles, JV Neel, The costs of human inbreeding and their implications for variations at the DNA level. Nat Genet 8, 117–121 (1994).
7
H Najmabadi, et al., Deep sequencing reveals 50 novel genes for recessive cognitive disorders. Nature 478, 57–63 (2011).
8
B Charlesworth, D Charlesworth, The genetic basis of inbreeding depression. Genet Res 74, 329–340 (1999).
9
NE Morton, JF Crow, HJ Muller, An estimate of the mutational damage in man from data on consanguineous marriages. Proc Natl Acad Sci USA 42, 855–863 (1956).
10
PK Joshi, et al., Directional dominance on stature and cognition in diverse human populations. Nature 523, 459–462 (2015).
11
M Kardos, G Luikart, FW Allendorf, Measuring individual inbreeding in the age of genomics: Marker-based measures are better than pedigrees. Heredity 115, 63–72 (2015).
12
MC Keller, et al., Runs of homozygosity implicate autozygosity as a schizophrenia risk factor. PLoS Genet 8, e1002656 (2012).
13
T Lencz, et al., Runs of homozygosity reveal highly penetrant recessive loci in schizophrenia. Proc Natl Acad Sci USA 104, 19942–19947 (2007).
14
MC Keller, PM Visscher, ME Goddard, Quantification of inbreeding due to distant ancestors and its detection using dense single nucleotide polymorphism data. Genetics 189, 237–249 (2011).
15
S Gazal, et al., Inbreeding coefficient estimation with dense SNP data: Comparison of strategies and application to HapMap III. Hum Hered 77, 49–62 (2014).
16
S Purcell, et al., PLINK: A tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81, 559–575 (2007).
17
J Yang, SH Lee, ME Goddard, PM Visscher, GCTA: A tool for genome-wide complex trait analysis. Am J Hum Genet 88, 76–82 (2011).
18
A Gusev, et al., Partitioning heritability of regulatory and cell-type-specific variants across 11 common diseases. Am J Hum Genet 95, 535–552 (2014).
19
J Yang, et al., Genetic variance estimation with imputed variants finds negligible missing heritability for human height and body mass index. Nat Genet 47, 1114–1120 (2015).
20
D Speed, G Hemani, MR Johnson, DJ Balding, Improved heritability estimation from genome-wide SNPs. Am J Hum Genet 91, 1011–1021 (2012).
21
SH Lee, et al., Estimation of SNP heritability from dense genotype data. Am J Hum Genet 93, 1151–1155 (2013).
22
BK Bulik-Sullivan, et al., LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat Genet 47, 291–295 (2015).
23
MR Robinson, et al., Population genetic differentiation of height and body mass index across Europe. Nat Genet 47, 1357–1362 (2015).
24
DP Howrigan, et al., Genome-wide autozygosity is associated with lower general cognitive ability. Mol Psychiatry 21, 837–843 (2016).
25
Z Zhu, et al., Dominance genetic variation contributes little to the missing heritability for human complex traits. Am J Hum Genet 96, 377–385 (2015).
26
JP Mackenbach, Health and deprivation. Inequality and the North. Health Policy 10, 207 (1988).
27
AR Wood, et al., Defining the role of common variation in the genomic and biological architecture of adult human height. Nat Genet 46, 1173–1186 (2014).
28
J Yang, et al., Common SNPs explain a large proportion of the heritability for human height. Nat Genet 42, 565–569 (2010).
29
A Robert, B Toupance, M Tremblay, E Heyer, Impact of inbreeding on fertility in a pre-industrial population. Eur J Hum Genet 17, 673–681 (2009).
30
AH Bittles, JC Grant, SG Sullivan, R Hussain, Does inbreeding lead to decreased human fertility? Ann Hum Biol 29, 111–130 (2002).
31
MA Woodley, Inbreeding depression and IQ in a study of 72 countries. Intelligence 37, 268–276 (2009).
32
M Lynch, The genetic interpretation of inbreeding depression and outbreeding depression. Evolution 45, 622–629 (1991).
33
Z Szpiech, et al., Long runs of homozygosity are enriched for deleterious variation. Am J Hum Genet 93, 90–102 (2013).
34
N Allen, et al., UK Biobank: Current status and what it means for epidemiology. Health Policy Technol 1, 123–126 (2012).
35
J O’Connell, et al., Haplotype estimation for biobank-scale data sets. Nat Genet 48, 817–820 (2016).
Information & Authors
Information
Published in
Classifications
Submission history
Published online: July 26, 2017
Published in issue: August 8, 2017
Keywords
Acknowledgments
This research was supported by the Australian Research Council (DP130102666 and DP160103860), the Australian National Health and Medical Research Council (1078037 and 1107258), the National Institute of Health (NIH Grants R01AG042568 and P01GM099568), and the Sylvia and Charles Viertel Charitable Foundation. This research has been conducted using the UK Biobank Resource under Project 12505.
Notes
This article is a PNAS Direct Submission.
Authors
Competing Interests
The authors declare no conflict of interest.
Metrics & Citations
Metrics
Citation statements
Altmetrics
Citations
Cite this article
114 (32) 8602-8607,
Export the article citation data by selecting a format from the list below and clicking Export.
Cited by
Loading...
View Options
View options
PDF format
Download this article as a PDF file
DOWNLOAD PDFLogin options
Check if you have access through your login credentials or your institution to get full access on this article.
Personal login Institutional LoginRecommend to a librarian
Recommend PNAS to a LibrarianPurchase options
Purchase this article to access the full text.