A non-zero variance of Tajima’s estimator for two sequences even for infinitely many unlinked loci

Abstract

The population-scaled mutation rate, $\theta$, is informative on the effective population size and is thus widely used in population genetics. We show that for two sequences and $n$ unlinked loci, the variance of Tajima’s estimator ($\hat{\theta }$), which is the average number of pairwise differences, does not vanish even as $n\to \infty$. The non-zero variance of $\hat{\theta }$ results from a (weak) correlation between coalescence times even at unlinked loci, which, in turn, is due to the underlying fixed pedigree shared by gene genealogies at all loci. We derive the correlation coefficient under a diploid, discrete-time, Wright–Fisher model, and we also derive a simple, closed-form lower bound. We also obtain empirical estimates of the correlation of coalescence times under demographic models inspired by large-scale human genealogies. While the effect we describe is small ($\mathrm{Var}\left[\hat{\theta }\right]∕{\theta }^2\approx \mathcal{O}\left(N_e^{-1}\right)$), it is important to recognize this feature of statistical population genetics, which runs counter to commonly held notions about unlinked loci.

Introduction

The population-scaled mutation rate, $\theta$, is defined as $4N_e\mu$, where $N_e$ is the effective population size and $\mu$ is the mutation rate per locus per generation (Wakeley, 2009). Two classic estimators were developed for $\theta$, Watterson’s (based on the number of segregating sites Watterson, 1975) and Tajima’s (based on the average number of pairwise differences Tajima (1983), Tajima (1989)). For a single pair of sequences, both estimators are identical (denoted here as $\hat{\theta }$) and equal to the number of differences between the sequences.

Increasing the number of sampled individuals has limited ability to improve these estimates of $\theta$, because shared ancestry reduces the number of independent branches on which mutations can arise (Rosenberg and Nordborg, 2002). Felsenstein (2006) showed that the variance of maximum likelihood estimates of $\theta$ decreases approximately logarithmically with the number of individuals sampled. In contrast, the variance decreases inversely with the number of independent loci. Thus, to increase the accuracy of estimates of $\theta$, it is generally more effective to increase the number of independent loci than the sample size at each locus (see also e.g., Edwards “bibausep Beerli (2000), Pluzhnikov “bibausep Donnelly (1996) and references within).

Consider a set of $n$ unlinked loci located on different (non-homologous) chromosomes. We show here that even as $n\to \infty$, the variance of the resulting estimate of $\theta$ does not converge to zero, in contrast to what we may have naïvely assumed. This behavior results from the fact that coalescence times, even at unlinked loci, are in fact weakly correlated, due to the sharing of the same fixed underlying pedigree across all loci (Wakeley et al., 2012). By conditioning on the number of shared genealogical common ancestors, we derive a simple approximate lower bound, as a function of $N_e$, on the variance of $\hat{\theta }$ (Sections 2 The relation of the variance of, 3 Modeling the effect of the shared pedigree).

Unlinked loci may also be sampled from the same chromosome, separated by an infinitely high recombination rate. The correlation of coalescence times in such a case is higher, as the two loci may travel together for the first few generations. Therefore, the extent of the correlation, and thereby, the variance of $\hat{\theta }$, also depend on the sampling configuration. In Section 4, we derive the correlation coefficient analytically, as a function of the configuration and the effective population size, using a diploid discrete time Wright–Fisher model (DDTWF). This model is an extension of the haploid DTWF model, previously advocated by Bhaskar et al. (2014) for the study of large samples from finite populations.

Our results for the variance of $\hat{\theta }$ were obtained under the Wright–Fisher demographic model. To shed light on the variance of $\hat{\theta }$ under more realistic demographic models, in Section 5 we run simulations based on real, large-scale human genealogical data (Kaplanis et al., 2017). The pedigrees inspired by different human populations differ from each other and from the Wright–Fisher pedigrees in a number of ways, for example in the variance of the relatedness of any two randomly chosen individuals. These differences lead to differences in the variance of $\hat{\theta }$ for each population, even if they have the same effective population size. Finally, we study some properties of linked sites in Section 6.

We note that the dependence of gene genealogies at unlinked loci has been previously recognized, most recently in the context of matching probabilities. Specifically, the probability of the genotypes of two individuals to match at two or more loci was computed under the Wright–Fisherand other models, and shown to differ from the product of the corresponding one-locus probabilities Laurie “bibausep Weir (2003), Song “bibausep Slatkin (2007), Bhaskar “bibausep Song (2009). In earlier literature, this effect was demonstrated in the context of identity-by-descent probabilities at unlinked loci (Weir and Cockerham, 1969) and implicitly in results on linkage disequilibrium (Ohta and Kimura, 1969). However, the treatment of this effect in the context of effective population size estimation is to our knowledge new.