Ratio estimators

The general problem of forming ratio estimators involves random variables a_i and b_i calculated from genotypes at each locus i, such that E[a_i] = Ac_i and E[b_i] = Bc_i and the goal is to estimate $$ \frac{A}{B}$$. A and B are constants shared across loci (given by F_ST or $$ {\phi }_{jk}^T$$), while c_i depends on the ancestral allele frequency $$ p_i^T$$ and varies per locus. The problem is that the single-locus estimator $$ \frac{a_i}{b_i}$$ is biased, since $$ E\left[\frac{a_i}{b_i}\right]\ne \frac{E\left[a_i\right]}{E\left[b_i\right]}=\frac{A}{B}$$, which applies to ratio estimators in general [69]. Below we study two estimator families that combine large numbers of loci to better estimate $$ \frac{A}{B}$$.

Convergence

The solution we recommend is the “ratio-of-means” estimator $$ \frac{{\widehat{A}}_m}{{\widehat{B}}_m}$$, where $$ {\widehat{A}}_m=\frac{1}{m}{\sum }_{i=1}^m{a}_i$$, and $$ {\widehat{B}}_m=\frac{1}{m}{\sum }_{i=1}^m{b}_i$$, which is common for F_ST estimators [4, 17, 19, 23, 70]. Note that $$ E\left[{\widehat{A}}_m\right]=A{\overline{c}}_m$$ and $$ E\left[{\widehat{B}}_m\right]=B{\overline{c}}_m$$, where $$ {\overline{c}}_m=\frac{1}{m}{\sum }_{i=1}^m{c}_i$$. We will assume bounded terms (|a_i|, |b_i| ≤ C for some finite C), a convergent $$ {\overline{c}}_m\to c$$, and Bc ≠ 0, which are satisfied by common estimators. Given independent loci, we prove almost sure convergence to the desired quantity (S1 Text),

In order to test if a given ratio-of-means estimator converges to its ratio of expectations as in Eq (10), the following three conditions can be tested. (i) The expected values of each term a_i, b_i must be calculated and shown to be of the form E[a_i] = Ac_i and E[b_i] = Bc_i for some A and B shared by all loci i and some c_i that may vary per locus i but must be shared by both E[a_i], E[b_i]. In the estimators we study, A and B are functions of IBD probabilities such as $$ {\phi }_{jk}^T$$ and F_ST, while c_i is a function of $$ p_i^T$$ only. (ii) The mean c_i must converge to a non-zero value for infinite loci. (iii) Both |a_i|, |b_i| ≤ C must be bounded for all i by some finite C (the estimators we study usually have C = 1 or C = 4). If these conditions are satisfied, then Eq (10) holds for independent loci and the A and B found in the first step. See the next section for an example application of this procedure to an F_ST estimator.

Another approach is the “mean-of-ratios” estimator $$ \frac{1}{m}{\sum }_{i=1}^m\frac{a_i}{b_i}$$, used often to estimate kinship coefficients [16, 27, 30–35] and F_ST [46]. If each $$ \frac{a_i}{b_i}$$ is biased, their average across loci will also be biased, even as m → ∞. However, if $$ E\left[\frac{a_i}{b_i}\right]\to \frac{A}{B}$$ for all loci i = 1, …, m as the number of individuals n → ∞, and $$ Var\left(\frac{a_i}{b_i}\right)$$ is bounded, then

The F_ST estimator for independent subpopulations and infinite subpopulation sample sizes

The WC, Weir-Hill, and Hudson method-of-moments estimators have small sample size corrections that remarkably make them consistent (as the number of independent loci m goes to infinity) for finite numbers of individuals. However, these small sample corrections also make the estimators unnecessarily cumbersome for our purposes (see Methods, section Previous F_ST estimators for the independent subpopulations model for complete formulas). In order to illustrate clearly how these estimators behave, both under the independent subpopulations model and for arbitrary structure, here we construct simplified versions that assume infinite sample sizes per subpopulation (Methods, section Previous F_ST estimators for the independent subpopulations model). This simplification corresponds to eliminating statistical sampling, leaving only genetic sampling to analyze [71]. Note that our simplified estimator nevertheless illustrates the general behavior of the WC, Weir-Hill, and Hudson estimators under arbitrary structure, and the results are equivalent to those we would obtain under finite sample sizes of individuals. While the Hudson F_ST estimator compares two subpopulations [23], based on that work we derive a generalized “HudsonK” estimator for more than two subpopulations in Methods, section Generalized HudsonK F_ST estimator. Note that HudsonK, first derived in [58], also equals the Weir-Goudet F_ST estimator for subpopulations [21] when loci are biallelic, which was derived independently using allele matching (S1 Text).

Under infinite subpopulation sample sizes, the allele frequencies at each locus and every subpopulation are known. Let j ∈ {1, …, n} index subpopulations rather than individuals and π_ij be the true allele frequency in subpopulation j at locus i. Note that π_ij are not estimated allele frequencies, but rather true subpopulation allele frequencies; this abstraction does not result in a practical estimation approach, but it greatly simplifies understanding of bias for subpopulations in a setting where there there is no statistical sampling. Although in this analysis of F_ST estimators the π_ij values are applied to subpopulations, for coherence with our previous work we shall call them “individual-specific allele frequencies” (IAF) [60, 61]. Whether for individuals or subpopulations, the key assumption is that IAFs satisfy the coancestry model of Eq (6). In this special case of infinite subpopulation sample sizes, all of WC, Weir-Hill, and HudsonK simplify to the following F_ST estimator for independent subpopulations (“indep”; derived in Methods, section Previous F_ST estimators for the independent subpopulations model):

F_ST estimation under the independent subpopulations model

Under the independent subpopulations model $$ {\theta }_{jk}^T=0$$ for j ≠ k, where T is the most recent common ancestor (MRCA) population of the set of subpopulations. Note that the estimator in Eq (13) can be derived directly from Eq (6) and these assumptions using the method of moments (ignoring the existence of previous F_ST estimators; S1 Text). The expectations of the two recurrent terms in Eq (13) are

Before applying the convergence result in Eq (10), we test that the three conditions listed in section Assessing the accuracy of genome-wide ratio estimators are met. Condition (i): The locus i terms are $$ a_i={\widehat{\sigma }}_i^2$$ and $$ b_i={\widehat{p}}_i^T(1-{\widehat{p}}_i^T)+\frac{1}{n}{\widehat{\sigma }}_i^2$$, which satisfy E[a_i] = Ac_i and E[b_i] = Bc_i with A = F_ST, B = 1, and $$ c_i=p_i^T(1-p_i^T)$$. Condition (ii): $$ {\overline{c}}_m\to c=E\left[p_i^T\right(1-p_i^T\left)|T\right]\ne 0$$ over the $$ p_i^T$$ distribution across loci. Condition (iii): Since $$ 0\le {\pi }_{ij},{\widehat{p}}_i^T\le 1$$, then $$ 0\le {\widehat{\sigma }}_i^2\le 1$$ and $$ 0\le {\widehat{p}}_i^T(1-{\widehat{p}}_i^T)\le \frac{1}{4}$$, and since n ≥ 2, C = 1 bounds both |a_i| and |b_i|. Therefore, for independent loci,

F_ST estimation under arbitrary coancestry

Now we consider applying the independent subpopulations F_ST estimator to dependent subpopulations. The key difference is that now $$ {\theta }_{jk}^T\ne 0$$ for every (j, k) will be assumed in our coancestry model in Eq (6), and now T may be either the MRCA population of all subpopulations or a more ancestral population. In this general setting, (j, k) may index either subpopulations or individuals. The two terms of $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ now satisfy

The F_ST estimator for independent subpopulations converges more generally to

Since $$ \frac{1}{n}{F}_{\text{ST}}\le {\overline{\theta }}^T\le F_{\text{ST}}$$ (S1 Text), this estimator has a downward bias in the general setting: it is asymptotically unbiased ($$ {\widehat{F}}_{\text{ST}}^{\text{indep}}\underset{m\to \infty }{\overset{\text{a.s.}}{\to }}{F}_{\text{ST}}$$) only when $$ {\overline{\theta }}^T=\frac{1}{n}{F}_{\text{ST}}$$, while bias is maximal when $$ {\overline{\theta }}^T=F_{\text{ST}}$$, where $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}\underset{m\to \infty }{\overset{\text{a.s.}}{\to }}0$$. For example, if $$ min{\theta }_{jk}^T=0$$ so the MRCA population T is fixed, but n is large and $$ {\theta }_{jk}^T\approx F_{\text{ST}}$$ for most pairs of subpopulations, then $$ {\overline{\theta }}^T\approx F_{\text{ST}}$$ as well, and $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}\approx 0$$. Therefore, the magnitude of the bias of $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ is unknown if $$ {\overline{\theta }}^T$$ is unknown, and small $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ estimates may arise even if F_ST is very large.

Coancestry estimation as a method of moments

Since the generalized F_ST is given by coancestry coefficients $$ {\theta }_{jj}^T$$ in Eq (9), a new F_ST estimator could be derived from estimates of $$ {\theta }_{jj}^T$$. Here we attempt to define a method-of-moments estimator for $$ {\theta }_{jk}^T$$, and find an underdetermined estimation problem, just as for F_ST. This is consistent with IBD parameters in general requiring a reference population to be determined [40], whereas in this subsection this reference population is unspecified.

Given IAFs and the coancestry model of Eq (6), the first and second moments that average across loci are

Suppose first that only $$ {\theta }_{jj}^T$$ are of interest. There are n estimators given by Eq (16) with j = k, each corresponding to an unknown $$ {\theta }_{jj}^T$$. However, all these estimators share two nuisance parameters: $$ {\overline{p}}^T$$ and $$ {\overline{p^2}}^T$$. While $$ {\overline{p}}^T$$ can be estimated from Eq (15), there are no more equations left to estimate $$ {\overline{p^2}}^T$$, so this system is underdetermined. The estimation problem remains underdetermined if all $$ \frac{n(n+1)}{2}$$ estimators in Eq (16) are considered rather than only the j = k cases. Therefore, we cannot estimate coancestry coefficients consistently using only the first two moments without additional assumptions.

Characterization of the standard kinship estimator

Here we analyze a standard kinship estimator that is frequently used [16, 27, 30–36]. We generalize this estimator to use weights in estimating the ancestral allele frequencies, and we write it as a ratio-of-means estimator due to the favorable theoretical properties of this format as detailed in the earlier section Assessing the accuracy of genome-wide ratio estimators:

Utilizing the kinship model for genotypes from Eq (5), we find that Eq (18) converges to

The estimators in Eqs (18) and (21) are unbiased when $$ {\widehat{p}}_i^T$$ is replaced by $$ p_i^T$$ [16, 32, 36], and are consistent when $$ {\widehat{p}}_i^T$$ is consistent [60]. Surprisingly, $$ {\widehat{p}}_i^T$$ in Eq (17) is not consistent (it does not converge almost surely to $$ p_i^T$$) for arbitrary population structures, which is at the root of the bias in Eqs (19) to (21). In particular, although $$ {\widehat{p}}_i^T$$ is unbiased, its variance (S1 Text, and some special cases shown elsewhere, e.g., [19]),

Estimation of coancestry coefficients from IAFs

Here we form a coancestry coefficient estimator analogous to Eq (18) but using IAFs. Assuming the moments in Eq (6), this estimator and its limit are

F_ST estimator based on the standard kinship estimator

Since the generalized F_ST is defined as a mean inbreeding coefficient in Eq (3), here we study the F_ST estimator constructed as $$ {\widehat{F}}_{\text{ST}}^{\text{std}}={\sum }_{j=1}^n{w}_j{\widehat{f}}_j^{T,\text{std}}$$ where $$ {\widehat{f}}_j^{T,\text{std}}$$ is the inbreeding estimator derived from the standard kinship estimator. Although $$ {\widehat{f}}_j^{T,\text{std}}$$ is biased, we nevertheless plug it into our definition of F_ST so that we may study how bias manifests. Note that we do not recommend utilizing this F_ST estimator in practice, but we find these results informative for identifying how to proceed in deriving new estimators in the following section.

Remarkably, the three $$ f_j^T$$ estimators in Eqs (20) and (21) give exactly the same plug-in $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ if the weights in F_ST and $$ {\widehat{p}}_i^T$$ in Eq (17) match, namely

Like $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ in Eq (13), $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eqs (25) and (26) are downwardly biased since $$ 0\le {\overline{\phi }}^T,{\overline{\theta }}^T$$. $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eq (26) may converge arbitrarily close to zero since $$ {\overline{\theta }}^T$$ can be arbitrarily close to F_ST (S1 Text). Moreover, although $$ {\overline{\phi }}^T\approx {\overline{\theta }}^T$$ for large n (see Eq (8) and panel “Coancestry in Terms of Kinship” in Fig 1), in extreme cases $$ {\overline{\phi }}^T$$ can exceed F_ST under the coancestry model (where $$ {\overline{\theta }}^T\le {\overline{\phi }}^T$$) and also under extreme local kinship, where $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eq (25) converges to a negative value.

Adjusted consistent oracle F_ST estimators and the “bias coefficient”

Here we explore two adjustments to $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ from IAFs in Eq (26) that rely on having minimal additional information needed to correct its bias. These “oracle” approaches require information that is not known in practice, but this exercise helps us understand the problem more deeply and finds further connections between the various F_ST estimators.

If $$ {\overline{\theta }}^T$$ is known, the bias in Eq (26) can be reversed, yielding the consistent estimator

Lastly, using either Eqs (26) or (29), the relative error of $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ converges to

General approach

In this subsection we develop our new estimator in two steps. First, we compute a new statistic A_jk that is proportional in the limit of infinite loci to $$ {\phi }_{jk}^T-1$$ times a nuisance factor v^T. Second, we estimate and remove v^T to yield the proposed estimator $$ {\widehat{\phi }}_{jk}^{T,\text{new}}$$. $$ {\widehat{A}}_{\text{min}}$$—an estimator of the limit of the minimum A_jk—yields v^T if the least related pair of individuals in the data has $$ {\phi }_{jk}^T=0$$, which sets T to the MRCA population of all the individuals in the data. The new kinship estimator immediately results in new inbreeding ($$ {\widehat{f}}_j^{T,\text{new}}$$) and F_ST ($$ {\widehat{F}}_{\text{ST}}^{\text{new}}$$) estimators. This general approach leaves the implementation of $$ {\widehat{A}}_{\text{min}}$$ open; the simple implementation applied in this work is described in subsection Proof-of-principle kinship estimator using subpopulation labels, but our method can be readily improved by substituting in a better $$ {\widehat{A}}_{\text{min}}$$ in the future.

Applying the method of moments to Eq (5), we derive the following statistic (S1 Text), whose expectation is proportional to $$ {\phi }_{jk}^T-1$$:

In general, suppose $$ {\widehat{A}}_{\text{min}}$$ is a consistent estimator of the limit of the minimum E[A_jk|T], or equivalently,

Proof-of-principle kinship estimator using subpopulation labels

To showcase the potential of the new estimators, we implement a simple proof-of-principle version of $$ {\widehat{A}}_{\text{min}}$$ needed for our new kinship estimator ($$ {\widehat{\phi }}_{jk}^{T,\text{new}}$$ in Eq (34)). This $$ {\widehat{A}}_{\text{min}}$$ relies on an appropriate partition of the n individuals into K subpopulations (denoted S_u for u ∈ {1, …, K}), where the only requirement is that the kinship coefficients between pairs of individuals across the two most unrelated subpopulations is zero, as detailed below. Note that, unlike the the independent subpopulations model of section F_ST estimation based on the independent subpopulations model, these K subpopulations need not be independent nor unstructured. The desired estimator $$ {\widehat{A}}_{\text{min}}$$ is the minimum average A_jk over all subpopulation pairs:

This estimator should work well for individuals truly divided into subpopulations, but may be biased for a poor choice of subpopulations, in particular if the minimum mean $$ {\phi }_{jk}^T$$ between subpopulations is far greater than zero. For this reason, inspection of the kinship estimates is required and careful construction of appropriate subpopulations may be needed. See our analysis of human data for detailed examples [59]. Future work could focus on a more general $$ {\widehat{A}}_{\text{min}}$$ that circumvents the need for subpopulations of our proof-of-principle estimator.

Comparison to the Weir-Goudet kinship estimator for individuals

Here we analyze the Weir-Goudet (WG) kinship estimator for individuals [20–22]. This has connections to our new estimator but differs in having the goal of estimating linearly-transformed kinship values. In our framework, the WG estimator is given by

Note that the original WG does not estimate F_ST from individuals as considered above; instead, F_ST is estimated from coancestry estimates for subpopulations (which equals our HudsonK for biallelic loci, S1 Text) [20–22]. For completeness, we consider both kinds of F_ST estimates in the evaluations that follow.

Overview of simulations

We simulate genotypes from two models to illustrate our results when the true population structure parameters are known. Both simulations have clearly-defined IBD probability parameters in terms of the MRCA population. The first simulation satisfies the independent subpopulations model that existing F_ST estimators assume. The second simulation is from an admixture model with no independent subpopulations and pervasive kinship designed to induce large downward biases in existing kinship and F_ST estimators (Fig 2). This admixture scenario resembles the population structure we estimated for Hispanics in the 1000 Genomes Project [59]: compare the simulated kinship matrix (Fig 2B) and admixture proportions (Fig 3C) to our estimates on the real data [59]. Both simulations have n = 1000 individuals, m = 300, 000 loci, and K = 10 subpopulations or intermediate subpopulations. These simulations have F_ST = 0.1, comparable to previous estimates between human populations (in 1000 Genomes, the estimated F_ST between CEU (European-Americans) and CHB (Chinese) is 0.106, between CEU and YRI (Yoruba from Nigeria) it is 0.139, and between CHB and YRI it is 0.161 [23]).

Fig 2

Coancestry matrices of simulations.

Both panels have n = 1000 individuals along both axes, K = 10 subpopulations (final or intermediate), and F_ST = 0.1. Color corresponds to $$ {\theta }_{jk}^T$$ between individuals j and k (equal to $$ {\phi }_{jk}^T$$ off-diagonal, $$ f_j^T$$ along the diagonal). (A) The independent subpopulations model has $$ {\theta }_{jk}^T=0$$ between subpopulations, and varying $$ {\theta }_{jj}^T$$ per subpopulation, resulting in a block-diagonal coancestry matrix. (B) Our admixture scenario models a 1D geography with extensive admixture and intermediate subpopulation differentiation that increases with distance, resulting in a smooth coancestry matrix with no independent subpopulations (no $$ {\theta }_{jk}^T=0$$ between blocks). Individuals are ordered along each axis by geographical position.

Fig 3

1D admixture scenario.

We model a 1D geography population that departs strongly from the independent subpopulations model. (A) K = 10 intermediate subpopulations, evenly spaced on a line, evolved independently in the past with F_ST increasing with distance, which models a sequence of increasing founder effects (from left to right) to mimic the global human population. (B) Once differentiated, individuals in these intermediate subpopulations spread by random walk modeled by Normal densities. (C) n = 1000 individuals, sampled evenly in the same geographical range, are admixed proportionally to the previous Normal densities. Thus, each individual draws most of its alleles from the closest intermediate subpopulation, and draws the fewest alleles from the most distant populations. Long-distance random walks of intermediate subpopulation individuals results in kinship for admixed individuals that decays smoothly with distance in Fig 2B. (D) For F_ST estimators that require a partition of individuals into subpopulations, individuals are clustered by geographical position (K = 10).

The independent subpopulations simulation satisfies the HudsonK and BayeScan estimator assumptions: each independent subpopulation S_u has a different F_ST value of $$ f_{S_u}^T$$ relative to the MRCA population T (Fig 2A). Ancestral allele frequencies $$ p_i^T$$ are drawn uniformly between 0.01 and 0.5. Allele frequencies $$ p_i^{S_u}$$ for S_u and locus i are drawn independently from the Balding-Nichols (BN) distribution [3] with parameters $$ p_i^T$$ and $$ f_{S_u}^T$$. Every individual j in subpopulation S_u draws alleles randomly with probability $$ p_i^{S_u}$$. Subpopulation sample sizes are drawn randomly (Methods, section Simulations).

The admixture simulation corresponds to a “BN-PSD” model [6, 27, 34, 60, 77]: the intermediate subpopulations are independent subpopulations that draw $$ p_i^{S_u}$$ from the BN model, then each individual j constructs its allele frequencies as $$ {\pi }_{ij}={\sum }_{u=1}^K{p}_i^{S_u}{q}_{ju}$$, which is a weighted average of the subpopulation allele frequencies $$ p_i^{S_u}$$ with the admixture proportions q_ju of individual j and subpopulation u as weights (which satisfy $$ {\sum }_{u=1}^K{q}_{ju}=1$$), as in the Pritchard-Stephens-Donnelly (PSD) admixture model [63–65]. We constructed q_ju that model admixture resulting from spread by random walk of the intermediate subpopulations along a one-dimensional geography, as follows. Intermediate subpopulations S_u are placed on a line with differentiation $$ f_{S_u}^T$$ that grows with distance, which corresponds to a serial founder effect (Fig 3A). Upon differentiation, individuals in each S_u spread by random walk, a process modeled by Normal densities (Fig 3B). Admixed individuals derive their ancestry proportional from these Normal densities, resulting in a genetic structure governed by geography (Figs 3C and 2B) and departing strongly from the independent subpopulations model (Fig 3D). The amount of spread—which sets the mean kinship across all individuals—was chosen to give a bias coefficient of $$ s^T=\frac{{\overline{\theta }}^T}{F_{\text{ST}}}=0.5$$, which by Eq (32) results in a large downward bias for $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ (in contrast, the independent subpopulations simulation has s^T = 0.1). The true coancestry and F_ST parameters of this simulation are given by the $$ f_{S_u}^T$$ values of the intermediate subpopulations and the admixture coefficients q_ju of the individuals via the following equations [57]:

Evaluation of F_ST estimators

Our admixture simulation illustrates the large biases that can arise if existing F_ST estimators that require independent subpopulations or F_ST estimates derived from existing kinship estimators are misapplied to arbitrary population structures to estimate the generalized F_ST, and demonstrate the higher accuracy of our new F_ST estimator ($$ {\widehat{F}}_{\text{ST}}^{\text{new}}$$ given by the combination of Eqs (36) and (37)). The WC F_IT (total inbreeding) estimator was also evaluated.

First, we test these estimators in our independent subpopulations simulation. The HudsonK (Methods, section Generalized HudsonK F_ST estimator) and BayeScan F_ST estimators are consistent in this simulation, since their assumptions are satisfied (Fig 4A). The WC F_ST estimator assumes that $$ f_{S_u}^T=F_{\text{ST}}$$ for all subpopulations S_u, which does not hold; nevertheless, WC has only a small bias (Fig 4A). The WC F_IT estimator arrives at similar estimates, as it should since there is no local inbreeding, so the true F_IT also equals F_ST. The Weir-Hill estimator permits different $$ f_{S_u}^T$$ values per subpopulation, but assigns equal weight to individuals rather than subpopulations (Methods, section The Weir-Hill F_ST estimator), resulting in a slightly different target F_ST (we verified that these estimates are unbiased for this F_ST). For comparison, we show the standard kinship-based $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eq (25) (weights from Methods, section Simulations) and $$ {\widehat{F}}_{\text{ST}}^{\text{WG}}$$ based on the Weir-Goudet kinship estimates for individuals, both of which do not have corrections that would make them consistent under the independent subpopulations model. Since the number of subpopulations K is large, $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ has a small relative bias of about $$ s^T=\frac{1}{K}=10\%$$ (Fig 4A); greater bias is expected for smaller K. Our new F_ST estimator has a very small bias in this simulation resulting from estimating the minimum kinship from the smallest kinship between subpopulations (see Eq (37)) rather than their average as HudsonK does implicitly (Fig 4A).

Fig 4

Evaluation of F_ST estimators.

The Weir-Cockerham, Weir-Hill, Weir-Goudet (for individuals), HudsonK (equal to Weir-Goudet for subpopulations, S1 Text), BayeScan, $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eq (25) derived from the standard kinship estimator, and our new F_ST estimator in Eqs (34) and (37), are evaluated on simulated genotypes from our two models (Fig 2). The Weir-Cockerham F_IT estimator was also included to show that estimation of total inbreeding behaves similarly to F_ST estimators. (A) The independent subpopulations model required by the Weir-Hill, HudsonK, and BayeScan F_ST estimators. All but standard kinship ($$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$) and Weir-Goudet (for individuals) recover the target F_ST IBD probability in Eq (9) (red line) with small errors. (B) Our admixture scenario, which has no independent subpopulations, was constructed so $$ {\widehat{F}}_{\text{ST}}^{\text{std}}\approx \frac{1}{2}{F}_{\text{ST}}$$. Only our new estimates are accurate. The rest of these estimators give values smaller than the target F_ST IBD probability, which result from treating kinship as zero between every subpopulations imposed by geographic clustering (or between individuals for Standard Kinship and Weir-Goudet). The $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ estimator limit in Eq (14) (green dotted line) overlaps the true F_ST (red line) in (A) but not (B). Estimates (blue) include 95% prediction intervals (often too narrow to see) from 39 independently-simulated genotype matrices for each model (Methods, section Prediction intervals).

Next we test these estimators in our admixture simulation. To apply the F_ST estimators that require subpopulations to the admixture model, individuals are clustered into subpopulations by their geographical position (Fig 3D). We find that estimates of all existing methods are smaller than the true F_ST by nearly half, as predicted by the limit of $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ in Eq (14) (Fig 4B). The WC F_IT estimator obtains slightly larger estimates than the WC F_ST estimator, but overall remains as biased as the other F_ST estimators, showing that the use of a total inbreeding estimator for independent subpopulations displays the same bias as the corresponding F_ST estimator. By construction, the kinship-based $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ also has a large relative bias of about s^T = 50%; remarkably, all existing F_ST estimators for subpopulations suffer from comparable biases. Thus, the corrections for independent subpopulations present in the WC, Weir-Hill and HudsonK estimators, or the Bayesian likelihood modeling of BayeScan, are insufficient for accurate estimation of the target generalized F_ST (Eq (9)) in this admixture scenario. Only our new F_ST estimator achieves accurate estimates of the generalized F_ST in the admixture simulation (Fig 4B).

Evaluation of kinship estimators

Our admixture simulation illustrates the distortions of the standard kinship estimator $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ in Eq (18), the linearly-transformed kinship values given by the Weir-Goudet estimator, and demonstrates the improved accuracy of our new kinship estimator $$ {\widehat{\phi }}_{jk}^{T,\text{new}}$$ given by the combination of Eqs (34) and (37). Kinship matrix estimates and their limits are visualized as heatmaps in Fig 5, whereas estimator accuracy is shown directly in Fig 6. The limit of the standard estimator $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ in Eq (18) would have had a uniform bias if $$ {\overline{\phi }}_j^T={\overline{\phi }}^T$$ held for all individuals j. For that reason, our admixture simulation has varying differentiation $$ f_{S_u}^T$$ per intermediate subpopulation S_u (Fig 3A), which causes large differences in $$ {\overline{\phi }}_j^T$$ per individual j and therefore large distortions in $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$. The Weir-Goudet approach estimates the linearly-transformed kinship values calculated in Eq (38).

Fig 5

Evaluation of kinship estimators.

Observed accuracy for two existing kinship coefficient estimators is illustrated in our admixture simulation and contrasted to the nearly unbiased estimates of our new estimator. Plots show n = 1000 individuals along both axes, and color corresponds to $$ {\phi }_{jk}^T$$ between individuals j ≠ k and to $$ f_j^T$$ along the diagonal ($$ f_j^T$$ is in the same scale as $$ {\phi }_{jk}^T$$ for j ≠ k; plotting $$ {\phi }_{jj}^T$$, which have a minimum value of $$ \frac{1}{2}$$, would result in a discontinuity in this figure). (A) True kinship matrix. (B) Estimated kinship using our new estimator in Eqs (34) and (37) from simulated genotypes recovers the true kinship matrix with high accuracy. (C) Theoretical limit of $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ in Eq (19) as the number of independent loci goes to infinity demonstrates the accuracy of our bias predictions under the kinship model. (D) Standard kinship estimates $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ given by Eq (18) from simulated genotypes are downwardly biased on average and distorted by pair-specific amounts. (E) Theoretical limit of the Weir-Goudet kinship estimator given by Eq (38). (F) Weir-Goudet kinship estimates from the same simulated genotypes agree with our calculated limit.

Fig 6

Accuracy of kinship estimators.

Here the estimated kinship values are directly compared to their true values, in the same admixture simulation data (n = 1000 individuals) shown in the previous figure. (A) Kinship between different individuals (excluding inbreeding). The new estimator has practically no bias in this evaluation (falls on the 1-1 dashed gray line). The standard estimator has a complex, non-linear bias that covers a large area of errors. (B) Inbreeding comparison, shows the bias of the standard estimate follows a different pattern for inbreeding compared to kinship between individuals. To better visualize and compare data across panels, a random subset of n points (out of the original n(n − 1)/2 unique individual pairs) were plotted in (A), matching the number of individuals (number of points in (B)).

Our new kinship estimator (Fig 5B) recovers the true kinship matrix of this complex population structure (Fig 5A), with an RMSE of 2.83% relative to the mean $$ {\phi }_{jk}^T$$ (Fig 6). In contrast, estimates using the standard estimator have a large overall downward bias (Fig 5C), resulting in an RMSE of 115.72% from the true $$ {\phi }_{jk}^T$$ relative to the mean $$ {\phi }_{jk}^T$$ (Fig 6). Additionally, estimates from $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ are very distorted, with an abundance of $$ {\widehat{\phi }}_{jk}^{T,\text{std}}<{\phi }_{jk}^T$$ cases—some of which are negative estimates (blue in Fig 5C)—but remarkably also cases with $$ {\widehat{\phi }}_{jk}^{T,\text{std}}>{\phi }_{jk}^T$$ (top left corner of Figs 5C and 6).

Now we compare the convergence of the ratio-of-means and mean-of-ratios versions of the standard kinship estimator to their biased limit we calculated in Eq (19) (Fig 5D). The ratio-of-means estimate $$ {\widehat{\phi }}_{jk}^{T,\text{std}}$$ (Fig 5C) has an RMSE of 2.14% from its limit relative to the mean $$ {\phi }_{jk}^T$$. In contrast, the mean-of-ratios estimates that are prevalent in the literature have a greater RMSE of 10.77% from the same limit in Eq (19). Thus, as expected from our theoretical results in section Assessing the accuracy of genome-wide ratio estimators, the ratio-of-means estimate is much closer to the desired limit than the mean-of-ratio estimate. The distortions are similar for the estimator that uses IAFs in Eq (24), with reduced RMSEs from its limit of 0.32% and 8.82% for the ratio-of-means and mean-of-ratios estimates, respectively.

Evaluation of oracle-adjusted F_ST estimators

Here we verify additional calculations for the bias of the standard kinship-based estimator $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ and the unbiased adjusted “oracle” F_ST estimators that require the true mean kinship $$ {\overline{\phi }}^T$$ or the bias coefficient s^T to be known. Note that $$ {\widehat{F}}_{\text{ST}}^{\text{new}}$$ in Eq (36) is related but not identical to these oracle estimators. We tested both IAF (Fig 7A) and genotype (Fig 7B) versions of these estimators. The unadjusted $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ in Eq (26) is severely biased (blue in Fig 7) by construction, and matches the calculated limit for IAFs and genotypes (green lines in Fig 7, which are close because $$ {\overline{\phi }}^T\approx {\overline{\theta }}^T$$). In contrast, the two consistent adjusted estimators $$ {\widehat{F}}_{\text{ST}}^{\prime }$$ and $$ {\widehat{F}}_{\text{ST}}^{\prime \prime }$$ in Eqs (27) and (31) estimate F_ST quite well (blue predictions overlap the true F_ST red line in Fig 7). However, $$ {\widehat{F}}_{\text{ST}}^{\prime }$$ and $$ {\widehat{F}}_{\text{ST}}^{\prime \prime }$$ are oracle methods, since they require parameters ($$ {\overline{\phi }}^T$$, $$ {\overline{\theta }}^T$$, s^T) that are not known in practice.

Fig 7

Evaluation of standard and adjusted F_ST estimators.

The convergence values we calculated for the standard kinship plug-in and adjusted F_ST estimators are validated using our admixture simulation. All adjusted estimators are unbiased but are “oracle” methods, since the mean kinship ($$ {\overline{\phi }}^T$$), mean coancestry ($$ {\overline{\theta }}^T$$), or bias coefficient ($$ s^T=\frac{{\overline{\theta }}^T}{F_{\text{ST}}}$$ for IAFs, replaced by $$ \frac{{\overline{\phi }}^T}{F_{\text{ST}}}$$ for genotypes) are usually unknown. (A) Estimation from individual-specific allele frequencies (IAFs): $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ is the standard coancestry plug-in estimator in Eq (26); $$ {\widehat{F}}_{\text{ST}}^{\prime }$$ “Adj. $$ {\overline{\theta }}^T$$” is in Eq (27); $$ {\widehat{F}}_{\text{ST}}^{\prime \prime }$$ “Adj. s” is in Eq (31). (B) For genotypes, $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ is given in Eq (25), and the adjusted estimators use $$ {\overline{\phi }}^T$$ rather than $$ {\overline{\theta }}^T$$. Lines: true F_ST (red line), limits of biased estimators $$ {\widehat{F}}_{\text{ST}}^{\text{std}}$$ (green lines, which differ slightly per panel). Estimates (blue) include 95% prediction intervals (too narrow to see) from 39 independently-simulated genotype matrices for our admixture model (Methods, section Prediction intervals).

Prediction intervals were computed from estimates over 39 independently-simulated IAF and genotype matrices (Methods, section Prediction intervals). Estimator limits are always contained in these intervals because the number of independent loci (m = 300, 000) is sufficiently large. Estimates that use genotypes have wider intervals than estimates from IAFs; however, IAFs are not known in practice, and use of estimated IAFs might increase noise. Genetic linkage, not present in our simulation, will also increase noise in real data.

The Weir-Cockerham F_ST estimator

The Weir-Cockerham (WC) F_ST estimator [17] estimates the coancestry parameter θ^T shared by each of the n independent subpopulation in consideration. Let $$ {\widehat{h}}_{ij}$$ denote the fraction of heterozygotes in subpopulation j for locus i. The ratio-of-means WC F_ST estimator and its limit for independent subpopulations ($$ {\theta }_{jk}^T=0$$ for j ≠ k) with equal differentiation ($$ {\theta }_{jj}^T={\theta }^T$$) is

Now we simplify this estimator as the sample size of every subpopulation becomes infinite. First set the sample size of every subpopulation n_j equal to their mean $$ \overline{n}$$, which implies C² = 0 and

The Weir-Hill F_ST estimator

Weir and Hill developed new estimators for subpopulation-specific F_ST values and considered the effects of non-independent subpopulations [4]. However, these estimators target linearly-transformed F_ST values, and recover the F_ST defined in Eq (9) only when subpopulations are independent [4], so we group them here with other estimators that strictly assume independent subpopulations. For simplicity, here we only consider the global F_ST estimator; the estimators of the coancestry matrix of the subpopulations was found to have the same overall linear transformation [4]. In the limit of infinite subpopulation sizes, this estimator also converges to the asymptotic F_ST estimator for independent subpopulations ($$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$) discussed in the main text.

The Weir-Hill (WH) F_ST estimator, simplified here for biallelic loci but extended to average over loci, and its limit, are given by

The Hudson F_ST estimator

The Hudson pairwise F_ST estimator [23] measures the differentiation of two subpopulations (j, k). The estimator and its limit for two independent subpopulations ($$ {\theta }_{jk}^T=0$$) is

Generalized HudsonK F_ST estimator

Here we derive the “HudsonK” estimator (first made available in [58]), which generalizes the Hudson pairwise F_ST estimator in Eq (40) to n independent subpopulations. This estimator also equals the recent Weir-Goudet F_ST estimator for subpopulations [21] (for biallelic loci; S1 Text). Note that for independent subpopulations, the F_ST of all the subpopulations equals the mean pairwise F_ST of every pair of subpopulations:

Like WC and Weir-Hill, $$ {\widehat{F}}_{\text{ST}}^{\text{HudsonK}}$$ simplifies to $$ {\widehat{F}}_{\text{ST}}^{\text{indep}}$$ in Eqs (11) to (13) in the limit of infinite sample sizes n_j → ∞, where $$ {\widehat{p}}_{ij}\to {\pi }_{ij}$$ for every (i, j).

Construction of subpopulation allele frequencies

We simulate K = 10 subpopulations S_u and m = 300, 000 independent loci. Every locus i draws $$ p_i^T~\text{Uniform}(0.01,0.5).$$ We set $$ f_{S_u}^T=\frac{u}{K}\tau ,$$ where τ ≤ 1 tunes F_ST. For the independent subpopulations model, $$ F_{\text{ST}}=\frac{1}{K}{\displaystyle {\sum }_{u=1}^K{f}_{S_u}^T}=\frac{\tau (K+1)}{2K},$$ so $$ \tau =\frac{2K{F}_{\text{ST}}}{K+1}$$ gives the desired F_ST (τ ≈ 0.18 for F_ST = 0.1). For the admixture model, τ is found numerically (τ ≈ 0.90 for F_ST = 0.1; see last subsection). Lastly, $$ p_i^{S_u}$$ values are drawn from the Balding-Nichols distribution,

Random subpopulation sizes

We randomly generate sample sizes r = (r_u) for K subpopulations and $$ {\sum }_{u=1}^K{r}_u=n=1000$$ individuals, as follows. First, draw x ∼ Dirichlet (1, …, 1) of length K and r = round(n x). While $$ {min}_u{r}_u<\frac{n}{3K}$$, draw a new r, to prevent small subpopulations (they do not occur in real data). Due to rounding, $$ {\sum }_{u=1}^K{r}_u$$ may not equal n as desired. Thus, while $$ \delta =n-{\sum }_{u=1}^K{r}_u\ne 0$$, a random u is updated to r_u ← r_u + sgn(δ), which brings δ closer to zero at every iteration. Weights for individuals j in S_u are $$ w_j=\frac{1}{K{r}_u}$$ so the generalized F_ST matches $$ F_{\text{ST}}=\frac{1}{K}{\sum }_{u=1}^K{f}_{S_u}^T$$ from the independent subpopulations model (see section The generalized F_ST for arbitrary population structures), which HudsonK estimates.

Admixture proportions from 1D geography

We construct q_ju from random-walk migrations along a one-dimensional geography. Let x_u be the coordinate of intermediate subpopulation u and y_j the coordinate of a modern individual j. We assume q_ju is proportional to f(|x_u − y_j|), or

Choosing σ and τ

Here we find values for σ (controls q_jk) and τ (scales $$ f_{S_u}^T$$) that give $$ s^T=\frac{1}{2}$$ and F_ST = 0.1 in the admixture model. In our simulation, $$ w_j=\frac{1}{n}$$ and $$ f_{S_u}^T=\frac{u}{K}\tau $$, so applying those parameters to Eq (39) gives $$ {\theta }_{jk}^T=\frac{\tau }{K}{\sum }_{u=1}^{K}u{q}_{ju}{q}_{ku}$$ and $$ F_{\text{ST}}=\frac{\tau }{nK}{\sum }_{j=1}^n{\sum }_{u=1}^{K}u{q}_{ju}^2$$. Therefore,