PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of molbiolevolLink to Publisher's site
 
Mol Biol Evol. 2009 March; 26(3): 501–512.
Published online 2008 November 6. doi:  10.1093/molbev/msn254
PMCID: PMC2915906

An Unbiased Estimator of Gene Diversity in Samples Containing Related Individuals

Abstract

Gene diversity is sometimes estimated from samples that contain inbred or related individuals. If inbred or related individuals are included in a sample, then the standard estimator for gene diversity produces a downward bias caused by an inflation of the variance of estimated allele frequencies. We develop an unbiased estimator for gene diversity that relies on kinship coefficients for pairs of individuals with known relationship and that reduces to the standard estimator when all individuals are noninbred and unrelated. Applying our estimator to data simulated based on allele frequencies observed for microsatellite loci in human populations, we find that the new estimator performs favorably compared with the standard estimator in terms of bias and similarly in terms of mean squared error. For human population-genetic data, we find that a close linear relationship previously seen between gene diversity and distance from East Africa is preserved when adjusting for the inclusion of close relatives.

Keywords: heterozygosity, identity by descent, kinship coefficient

Introduction

Gene diversity, or expected heterozygosity, is a frequently used measure of genetic variation applied in diverse areas of population genetics. Together with its counterpart, gene identity or expected homozygosity, it has been used to quantify genetic variation in populations (Driscoll et al. 2002; Hoelzel et al. 2002), evaluate genetic divergence and population relationships (Nei 1973; Ramachandran et al. 2005), detect inbreeding (Li and Horvitz 1953), measure linkage disequilibrium (Ohta 1980; Sabatti and Risch 2002), and test for the influence of natural selection (Watterson 1978; Depaulis and Veuille 1998; Sabeti et al. 2002).

Consider a polymorphic locus with I distinct alleles and a population with parametric allele frequencies p1, p2, … ,pI, where pi [set membership] [0, 1] and ∑i = 1Ipi = 1. The term “gene diversity,” which is defined as

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx1_ht.jpg
(1)

was proposed by Nei (1973), though the use of equation (1) as a measure of diversity dates to considerably earlier (Gini 1912; Simpson 1949; Gibbs and Martin 1962).

Now consider a sample of n observations of alleles, in which the number of observations of allelic type i is ni. The count estimate of pi is p^i=ni/n. If no inbred or related individuals are included in the sample, then an unbiased estimator of gene diversity is (Nei and Roychoudhury 1974)

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx2_ht.jpg
(2)

If relatives or inbred individuals are included in the sample, then H^ is no longer an unbiased estimator of H. To understand why this statement is true, suppose that a sample contains a pair of close relatives. Because these individuals are related, they may share one or two alleles identically by descent (IBD) at a locus (compared with zero alleles shared IBD in unrelated individuals). As a result, estimation of pi is based on fewer independent observations than for a sample not containing any relatives. Although E[p^i]=pi when relatives are included, Var[p^i] is greater than it would be had no relatives been included. Observe that the computation of E[H^] involves a negative coefficient for E[p^i2]. Because E[p^i2]=Var[p^i]+E[p^i]2, E[H^] decreases as Var[p^i] increases. Thus, the inclusion of relatives results in a downward bias, so that E[H^]<H. For the case in which inbred unrelated individuals with known inbreeding coefficients are included in a sample, Weir (1989, 1996) provided the expectation of 1i=1Ip^i2, producing an unbiased estimator of gene diversity

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx3_ht.jpg
(3)

where An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx4_ht.jpg is the average inbreeding coefficient across individuals (see also Shete 2003). When inbred individuals are included, An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx4_ht.jpg≠0, and it follows that E[H^]<E[H^Weir]=H.

In this article, we conduct a detailed investigation of the case in which a sample includes related individuals. We derive an unbiased estimator of H for samples containing related individuals with known levels of relationship. Our derivation makes use of a formula of Bourgain et al. (2003) and McPeek et al. (2004) for the variance of count estimates of allele frequencies in samples containing inbred and related individuals. The resulting estimator incorporates kinship coefficients, the same quantitative descriptors of pairwise relationships that have been used in diverse problems involving relatives—such as evaluation of phenotypic covariances in families (Lange 2002), estimation of relatedness parameters (Weir et al. 2006), and quantitative-trait linkage analysis (Almasy and Blangero 1998). When a sample consists only of unrelated noninbred individuals, our new estimator H~ reduces to the standard estimator H^, and it reduces to H^Weir if inbred but not related individuals are included. Using data simulated based on allele frequencies from human populations, we find that the new estimator H~ corrects for bias generated by inclusion of related individuals and that it attains a mean squared error (MSE) comparable with that of H^. We apply this new estimator to microsatellite data from human population samples containing relatives and show that, compared with the standard estimator, it produces estimates closer to those obtained when excluding relatives from the analysis.

Theory

We assume that gene diversity is estimated from n/2 diploid individuals. Our aim is to obtain a bias-correction factor that can be incorporated into a new estimator of gene diversity, H~. We begin by computing Var[p^i] in a sample that may include relatives or inbred individuals. Var[p^i] was reported by Bourgain et al. (2003) and McPeek et al. (2004); we provide an alternative derivation that uses a generalization of the simpler method of Broman (2001). This approach was originally applied in a setting that did not consider inbreeding, and we generalize the computation to include inbreeding. Note that the variances of other estimators of allele frequencies have previously been derived in fairly general settings (McPeek et al. 2004) and that the estimator p^i is not a maximum likelihood estimator when related individuals are included in a sample (Boehnke 1991). However, our interest here is specifically on the count-based estimator of allele frequencies, as it is this estimator that is used in the standard estimator of gene diversity in equation (2).

Define Xk to be the number of alleles of type i that are carried by individual k at a particular locus. Xk can equal 0, 1, or 2, and E[Xk] = 2pi. Regardless of the relationships among individuals 1, 2, …, n/2, an unbiased estimator for pi, the frequency of allele i, is

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx5_ht.jpg
(4)

The variance of p^i is given by

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx6_ht.jpg
(5)

Suppose that individuals j and k are related. The coefficient of kinship between individuals j and k, Φj,k, is the probability that two alleles chosen at the locus—one from individual j and the other from individual k—are identical by descent. In the special case of j = k, the kinship coefficient is Φk,k = (1/2)(1 + fk), where fk is the inbreeding coefficient for individual k (Lange 2002, p. 81).

Conditional on the nature of the relationship between individuals j and k and on their inbreeding coefficients, the four alleles in the two individuals can take on one of nine condensed identity states (Jacquard 1974, p. 107). Let Δs = Pr[S = s], where the condensed identity state S ranges from 1 to 9 and the probability is conditional on the type of the relationship. Using table 1 and the fact that the kinship coefficient for the pair of individuals equals Δ1 + (1/2)(Δ3 + Δ5 + Δ7) + (1/4)Δ8 (Jacquard, 1974, p. 108), we obtain

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx7_ht.jpg

Because E[Xj] = E[Xk] = 2pi, it follows that

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx8_ht.jpg
(6)

Inserting the covariance into equation (5) yields

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx9_ht.jpg
(7)

where An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx10_ht.jpg is the average kinship coefficient across pairs of individuals (including comparisons of individuals with themselves). This result can be seen to be equivalent to the variance reported by McPeek et al. (2004, p. 361).

Proposition 1

Consider a locus with I distinct alleles, allele frequencies pi [set membership] [0, 1] and i=1Ipi=1. Suppose a sample from a population has n/2 possibly related and inbred individuals. Then an unbiased estimator for gene diversity is

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx11_ht.jpg
(8)

where Φj,k is the kinship coefficient of individuals j and k and An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx10_ht.jpg is the average kinship coefficient across pairs of individuals.

Proof

We need to show that E[H~]=H. Observing that E[p^i2]=Var[p^i]+E[p^i]2 and E[p^i]=pi, we apply equation (4) and then the variance of p^i in equation (7) to get

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx12_ht.jpg
Corollary 2

Consider a locus with I distinct alleles, allele frequencies pi [set membership] [0, 1] and i=1Ipi=1. Suppose a sample from a population has n/2 possibly related and inbred individuals. Let R be the set of distinct types of relative pairs in the sample. Further, let nR be the number of pairs of individuals with relationship type R and let ΦR be the kinship coefficient for each of these pairs. Then an unbiased estimator for gene diversity is

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx13_ht.jpg
(9)

where An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx14_ht.jpg is the average inbreeding coefficient across individuals and fk is the inbreeding coefficient for individual k.

Proof

Applying the definitions of An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx15_ht.jpg and Φk,k and the fact that Φj,k = 0 for a pair of “unrelated” individuals,

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx16_ht.jpg

Inserting this value for An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx15_ht.jpg into equation (8), we obtain the desired result.□

Table 1
Joint Distribution of the Numbers of i Alleles Carried by Individuals j and k Given Their Descent Configuration S, Assuming Allele i Has Frequency p

Note that if no related individuals are included in the sample, then R is the empty set, thus reducing H~ to H^Weir; if additionally no inbred individuals are included, then An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx4_ht.jpg=0 and H~ reduces to H^.

Corollary 3

Consider a locus with I distinct alleles, allele frequencies pi [set membership] [0, 1] and i=1Ipi=1. Suppose a sample from a population has n/2 noninbred individuals, among which q parent–offspring pairs, r full-sib pairs, and s second-degree (avuncular, grandparent–grandchild, and half-sib) relative pairs are included. Assuming the sample has no other relative pairs, an unbiased estimator for gene diversity is

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx17_ht.jpg
(10)
Proof

The kinship coefficients are ΦP = 1/4 for parent–offspring pairs, ΦF = 1/4 for full-sib pairs, and ΦS = 1/8 for second-degree pairs. If an individual k is not inbred, then fk = 0. For a sample without inbred individuals, An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx4_ht.jpg=0. Inserting the quantity and kinship coefficient for each of the three types of relative pairs into equation (9), we obtain equation (10).□

Corollary 4

Consider a locus with I distinct alleles, allele frequencies pi [set membership] [0, 1] and i=1Ipi=1. Suppose a sample from a population has n/2 possibly related and inbred individuals. Let R be the set of distinct types of relative pairs in the sample. Further, let nR be the number of pairs of individuals with relationship type R[set membership]R and let ΦR be the kinship coefficient for each of these pairs. Then the bias of H^ is always negative, increases in magnitude as H increases, and is given by

An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx18_ht.jpg
(11)

where An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx14_ht.jpg is the average inbreeding coefficient across individuals and fk is the inbreeding coefficient for individual k.

Proof

As shown in Corollary 2, H~=cH^, where c=n(n−1)/[n(n−1−An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx4_ht.jpg)−8∑R[set membership]RnRΦR]. Rearranging and taking the expected value gives E[H^]=E[H~]/c=H/c. The desired result follows from simplifying the expression for bias(H^), or (1 − c)H/c.□

Data from Human Populations

To examine the behavior of H~ in a realistic setting, we performed simulations and data analysis using microsatellite loci from the H1048 and H952 subsets (Rosenberg 2006) of the Human Genome Diversity Project–Centre d'Etude du Polymorphisme Humain (HGDP–CEPH) Cell Line Panel (Cann et al. 2002; Cavalli-Sforza 2005). The H1048 subset consists of 1,048 individuals in 53 populations. Among the 53 populations, the samples from 26 of them contain at least one pair of closely related individuals with either a first-degree (parent–offspring, full-sib) or second-degree (avuncular, grandparent–grandchild, and half-sib) relationship (table 2). The H952 subset is a collection of 952 individuals included in the larger H1048 subset. No two of the 952 individuals are believed to have a first- or second-degree relationship. Levels of relationship in H1048, as estimated previously from microsatellite genotypes (Rosenberg 2006), were treated here as known with certainty. Because no cycles were observed in pedigrees from the HGDP–CEPH panel (Rosenberg 2006), we assumed that none of the panel members were inbred. Genotypes at 783 autosomal microsatellite loci (Ramachandran et al. 2005; Rosenberg et al. 2005) were investigated in the H1048 and H952 data sets.

Table 2
The 26 Populations Containing Relatives in the H1048 Data Set (Modified from Rosenberg 2006, Supplementary tables 16 and 19)

Simulations

Simulation Procedure

Simulations based on the microsatellite loci were used to examine the properties of H~ and H^. For each of the 783 loci, we treated allele frequencies estimated from the H952 subset of individuals as true allele frequencies. The parametric gene diversity H was obtained for a locus as one minus the sum of the squares of these allele frequencies. All of our simulations assumed no inbreeding.

For a given locus, individual genotypes were simulated by sampling two alleles independently from the allele frequency distribution. To simulate a related individual with a given level of relationship to another individual, the number of alleles shared IBD with its relative was drawn under the appropriate probability distribution for the specified type of relative pair (parent–offspring, full-sib, or second-degree). This number of shared alleles (0, 1, or 2) was copied from a random individual that had already been generated and that had not yet been paired with a relative; if the number of alleles copied was 1, then an allele was chosen at random from the previously generated individual. The rest of the alleles, if any, were sampled independently from the allele frequency distribution. Gene diversity was estimated using H~ and H^ for samples with and without related individuals. We applied H^ both to entire samples as well to samples in which the “second” member of each relative pair was discarded. For each locus, simulated sets of individuals were obtained 100,000 times, and H^, H~, H^2, and H~2 were averaged across all replicates. The true value for gene diversity, H, was then subtracted from the mean of H^ and H~ to calculate bias for each estimator (and the result was squared to give bias squared). Variance of H^ was calculated by subtracting the square of the mean of H~ from the mean of H^2 (variance of H~ was calculated analogously). MSE was then calculated by adding variance and bias squared. Note that in our simulations, relative pairs were all disjoint, so that no individual was contained in multiple relative pairs; however, in our derivations, it is not required for relative pairs to be disjoint for H~ to be unbiased.

Simulation Results

To illustrate the performance of the estimators across the span of gene diversities present in the human microsatellite data set, loci were placed in increasing order by assumed parametric gene diversity, and six equally spaced loci—with the 112th, 224th, 336th, 448th, 560th, and 672nd highest values of gene diversity—were chosen for analysis. Similar results were obtained with all six loci (data not shown), and therefore, among the six loci only the locus with the lowest gene diversity (AAT263P, H = 0.6778) and the locus with the highest gene diversity (ACT3F12, H = 0.8263) were chosen for display. For both loci, table 3 shows the simulated MSE, variance, and bias squared for the different estimators, considering three different sample sizes and three combinations of the number of related individuals for each sample size. Because the simulation results are based on 100,000 replicate data sets, each of the quantities presented is small. However, it is possible to observe differences in the properties of the three estimators. Among the three estimators, H^ applied to full samples gives the lowest variance, H~ produces slightly higher variance, and H^ applied to samples with related individuals removed produces the highest variance. Bias squared is very close to zero for H^ applied to samples with related individuals removed, as well as for H~, but it is noticeably higher for H^ applied to full samples containing relatives. For the locus with the lower value of H (0.6778), H^ applied to full samples has the smallest MSE in all cases tested, although H~ has MSE very close to that of H^. However, for the locus with the higher value of H (0.8263), MSE is always smallest for H~. Therefore, H~ is not only unbiased, but it also has MSE comparable with—and sometimes smaller than—that of H^.

Table 3
MSE, Variance, and Bias Squared of Estimates for Data Simulated Based on Allele Frequencies at Two Loci (AAT263P and ACT3F12)

It is instructive to investigate the influence of specific variables on the MSE, variance, and bias squared of H~ and H^, by varying the simulation parameters over the space of gene diversities, sample sizes, and possible sets of relative pairs, and calculating MSE, variance, and bias squared for each scenario. We use H^full and H~full to denote H^ and H~ applied to a sample of individuals. For H^ applied to a sample in which related individuals are removed, we use the notation H^reduced.

Figure 1 displays the effect of sample size on MSE for each of the estimators, for scenarios in which all simulated individuals belong to relative pairs of a particular type. Here, the full and reduced samples consist of m and m/2 individuals, respectively. When q = m/2, r = m/2, or s = m/2, MSE is consistently lower for H^full and H~full (which have virtually identical MSE and therefore have overlapping lines in the graph) than for H^reduced. As the sample size increases, the MSEs of all estimators approach zero.

FIG. 1.
MSE as a function of sample size m for three different estimators. Each plot in a given row represents samples with a different type of relative pair. The numbers of parent–offspring, full-sib, and second-degree pairs are denoted by q, r, and ...

We next examined how the three estimators performed in simulated samples containing the same sample size and total number of relative pairs but with different combinations involving different numbers of parent–offspring, full-sib, and second-degree pairs. The same two loci that were analyzed in table 3 and figure 1 were investigated to show the effect of the combination of relative pairs at differing degrees of gene diversity. Figures 2 and and33 illustrate MSE, variance, and bias squared for each estimator as functions of the combination of types of relative pairs in a full sample of size 40 and a reduced sample of size 20 individuals. Each point in a triangle represents the number of parent–offspring, full-sib, and second-degree relative pairs in a sample; the sum of these quantities is equal to half the sample size. MSE and variance are always lower for H^full and H~full than for H^reduced, which relies on a smaller sample size, and H^full and H~full show similar trends. Bias squared for the unbiased H~full is similar to that for H^reduced, which eliminates relatives from the sample, whereas it is much larger for H^full. As the number of first-degree pairs is increased (decreasing the number of second-degree pairs), both variance and MSE increase. For H^full, as can be predicted from equation (11), bias squared also increases with an increase in the number of first-degree pairs. Because they are both unbiased estimators, H^full and H^reduced display no particular pattern for bias squared.

FIG. 2.
Heat maps of simulated MSE, variance, and bias squared for each estimator applied to a full sample of 40 and a reduced sample of 20 individuals, as functions of the mixture of types of relative pairs included in the sample. The simulation was based on ...
FIG. 3.
Heat maps of simulated MSE, variance, and bias squared for each estimator applied to a full sample of 40 and a reduced sample of 20 individuals, as functions of the mixture of types of relative pairs included in the sample. The simulation was based on ...

Finally, we studied the trends in MSE, variance, and bias squared for the estimators over the space of gene diversities, holding the full sample size fixed at 30 individuals and the reduced sample size fixed at 15. Unlike the analyses in table 3 and figures 13, which show results based on two representative loci, this analysis used simulations based on all 783 microsatellites. We considered a scenario in which the sample of 30 individuals consisted of 15 parent–offspring pairs. Figure 4 illustrates that for all three estimators, MSE and variance tend to decrease as gene diversity increases. Because H~full and H^reduced are both unbiased, bias squared shows no trend for these estimators. However, because bias for H^full is linear with respect to gene diversity (eq. 11), bias squared is quadratic. On the basis of equation (11), we predict An external file that holds a picture, illustration, etc.
Object name is molbiolevolmsn254fx19_ht.jpg, and a close match to this prediction was observed. The regression displayed in figure 4 has regression model [bias(H^full)]2=(7.187×105)H2.

FIG. 4.
MSE, variance, and bias squared for each estimator applied to a full sample of 30 and a reduced sample of 15 individuals, as functions of parametric gene diversity, considering simulated values based on each of the 783 loci. The simulations incorporated ...

Three main results can be observed in our simulations. First, H~ is unbiased and has comparable bias in samples containing relatives to that obtained by applying H^ to samples with relatives removed. Using H~, or excluding relatives and using H^, reduces the bias compared with using H^ without excluding relatives. Second, H^ has comparable (but consistently slightly higher) variance to the values obtained with H^ in samples containing relatives. Both H~ and H^ have lower variance in full samples of individuals than that of H^ in reduced samples that exclude relatives. Third, because H~ has less bias than H^ in samples containing relatives, H~ has comparable, and sometimes smaller, MSE to H^ (although its variance is larger). Both estimators have lower MSE than H^ applied to subsets that exclude relatives.

The properties of the estimators depend on a number of parameters. All estimators have lower MSE as sample size increases. In addition, the MSEs of H^ and H~ are smaller when second-degree relative pairs are investigated, in comparison to scenarios that include an equivalent number of first-degree pairs. Furthermore, the MSEs of H^ and H~ are generally smaller for loci with larger gene diversities, with the magnitude of the bias of H^ increasing linearly with increasing gene diversity.

We can conclude that for samples containing relatives, H~ has comparable variance to H^, with a considerable reduction of bias. H~ has comparable bias in a full sample to that of H^ applied to a reduced sample excluding relatives, with a considerable reduction of variance. Thus, H~ combines into a single estimator the desirable properties possessed by H^ applied to samples with relatives and by H^ applied to samples without relatives.

Application to Data

Notation

For convenience, we use the following notation: H^952 and H^1048 for application of H^ to the samples of 952 and 1,048 individuals, respectively, and H~952 and H~1048 for application of H~ to these samples. Note that because the H952 data set contains no relative pairs, H~952=H^952, and there is no need to consider H~952 separately. We also use the notation H^507, H^603, and H~603 when restricting our analysis to the 26 populations containing at least one relative pair; for each of the 27 remaining populations, the estimators H^ and H~ produce identical values.

Mean of the Estimator

For investigating the properties of H^ and H~ applied to the H1048 data set, because the true value of H is unknown for the actual data, we treated the value of H^952 for each locus as a substitute “true” value. Because H^ is unbiased when applied to data not containing relatives, H^952 provides a sensible proxy for the unknown true gene diversity. This approach enabled us to consider how estimates of H from data including relatives might differ from estimates based on the same data excluding all relatives. For each of the 53 populations, we computed the means of H^952, H^1048, and H~1048 across the 783 microsatellite loci. Because the true H is unknown and bias cannot be calculated, we instead examine the mean of H^1048 across loci minus the mean of H^952 across loci and the mean of H~1048 across loci minus the mean of H^952 across loci.

Figure 5 shows comparisons of the mean of H^1048H^952 across loci and the mean of H~1048H^952 across loci. In general, the three estimators produce similar estimates in a given population. However, notice that in figure 5A, H^1048 is reduced compared with H^952, a likely consequence of the bias of H^ when applied to samples containing relatives. When H~1048 is used in place of H^1048, because H~1048 corrects for the inclusion of known related individuals, there is a considerable reduction in the magnitude of the difference between the mean of the estimator (H^1048 or H~1048) across loci and the mean of H^952 across loci (fig. 5B). These observations are reflected in Wilcoxon signed rank tests that compare paired lists of mean heterozygosities across loci for the 53 populations (table 4). The P value for a comparison of H^1048 with H^952 was 8.804 × 10−6, suggesting that inclusion of relatives in a sample has a statistically significant impact on H^. In contrast, H~1048 and H^952 showed no significant difference, with a P value of 0.703 for the Wilcoxon signed rank test. Similar results were obtained for other comparisons of the three estimators. The mean across populations of H^952H~1048(3.262×104) was smaller than for H^952H^1048(2.387×103); the same was true for the mean of |H^952H~1048|(6.660×104) compared with the mean of |H^952H^1048|(2.387×103).

Table 4
Statistical Tests Applied to the Mean Gene Diversity across Loci
FIG. 5.
Comparison of the mean of H^1048H^952 and the mean of H~1048H^952. Each population is represented by a point colored based on the geographic location of the population, and the dotted line represents zero difference between the full-data ...

Comparable results were obtained when using only the 26 populations that contained relative pairs. The Wilcoxon signed rank test produced a statistically significant P value of 2.980 × 10−8 for H^603 compared with H^507 and a nonsignificant P value of 0.708 when comparing H~603 with H^507. The mean across populations of H^507H~603(6.649×104) was smaller than for H^507H^603(4.866×103), as was the mean of |H^507H~603|(1.358×103) relative to that of |H^507H^603|(4.866×103). In addition, similar numbers of populations had H~603>H^507 (12) and H~603<H^507 (14); by contrast, there were no populations with H^603>H^507.

Because estimators often have a trade-off between bias and variance, we investigated the relationship between the mean values across loci of H^603H^507 and H~603H^507 and the standard deviations of H^603 and H~603 across loci. We observed that compared with H^603, H~603 produces a noticeable decrease in the mean difference from H^507 with only a slight increase in the standard deviation (fig. 6). This result is somewhat analogous to the simulation-based result that H~ has less bias than H^ and comparable variance.

FIG. 6.
Comparison of the mean difference of an estimator (H^603 or H~603) from H^507 with the standard deviation of the estimator. Each population is represented by a point colored based on the geographic location of the population. Open and filled circles represent ...

Gene Diversity versus Distance from Africa

Based on an observed decline of gene diversity estimates with geographic distance from East Africa, Ramachandran et al. (2005) argued that the geographic expansion of modern humans can be described by a series of founder events originating in Africa. This analysis utilized the H^ estimator applied to the 783 microsatellites typed in the H1048 subset of individuals, excluding the Surui population. To evaluate how the results of Ramachandran et al. (2005) were affected by the bias of H^ in samples with close relatives, we analyzed the relationships of the three estimators of gene diversity—H^952, H^1048, and H~1048—with geographic distance from East Africa (fig. 7). Distance from Addis Ababa was measured in kilometers via waypoint routes and was based on the values from Rosenberg et al. (2005).

FIG. 7.
Gene diversity versus geographic distance from Addis Ababa, Ethiopia. (A) H^952 versus distance from Addis Ababa. The linear regression is given by H = 0.7778 − (7.955 × 10−6) × distance, with R2 = 0.856. (B) H^1048 versus ...

The three estimators produced relatively similar regressions (fig. 7), demonstrating that the close linear relationship of gene diversity and distance from Africa is not greatly affected by inclusion of relatives in the analysis. We observed very similar values for the coefficients of determination (R2) of linear regressions when using H^952, H^1048, and H~1048 (note that all three R2 values are higher than that reported by Ramachandran et al. (2005), whose lower value resulted from an error in the calculation of their fig. 4A). The Surui population, which has the smallest gene diversity and is the farthest population from Addis Ababa, deviates considerably from the regression line when using H^1048 to measure gene diversity (fig. 7B). When excluding the large number of relatives present in the Surui sample (H^952) or correcting for their inclusion (H~1048), the Surui population is not as extreme an outlier (fig. 7A and C).

Discussion

In this article, we have developed an unbiased estimator H~ for gene diversity in samples containing related and inbred individuals. The bias-correction factor in this estimator, which we derived from the variance of allele frequency estimates, depends only on the average kinship coefficient between pairs of sampled individuals. Using data simulated based on allele frequency distributions from human populations, we found that H~ performs well with regard to both bias and MSE. The bias generated by H~ applied to data including relatives is approximately the same as the bias generated by the standard estimator H^ applied to data containing only unrelated individuals. The MSE for H~ is comparable to—and often smaller than—the MSE of H^ when related individuals are included. Calculation of H~ relies only on sample allele frequencies and on the average kinship coefficient and is therefore easy to perform when relationships among individuals are known. Thus, the new estimator H~ offers a combination of unbiasedness, low MSE, and ease of computation, providing an improved approach to the estimation of gene diversity in samples containing relatives.

Using data from human populations, we found that H~ largely corrected a reduction in the standard estimator H^, producing estimates that were not significantly different from those obtained if we instead removed relatives from the data set and applied H^. This shift toward the values obtained in data without relatives occurred together with only a slight increase in standard deviation across loci relative to H^. However, by treating dependent observations as independent, H^ perhaps produces a smaller variance than is appropriate in samples with relatives. Thus, we conclude that as an alternative to removing relatives from samples containing relative pairs, H~ can be applied to obtain suitable gene diversity estimates.

When we applied H~ to the human data, a few populations still produced a “bias,” in that H~1048 remained considerably lower than H^952. The most noticeable of these populations are the Surui, Karitiana, and Pima populations from the Americas (fig. 5B); the “bias” was larger for these low-diversity populations, whereas theory predicts less bias when diversity is lower (eq. 11). It should first be noted that unlike for the other populations, inferences about second-degree relationships obtained by Rosenberg (2006) were somewhat uncertain for the Surui and Karitiana populations. Thus, table 2 and our analysis did not include inferred second-degree relationships in those populations, when in fact many are likely to be present. This is a likely reason why the “bias” in the Surui and Karitiana populations was only partially eliminated. For the Pima population, a likely explanation is that the sample contains many related individuals in extended families (Rosenberg, 2006), and our computation only adjusted for first- and second-degree relative pairs. If these higher order relationships had been fully known, however, it would have been possible to use our estimator to adjust for them.

Our estimator adjusts for inbreeding by averaging over inbreeding coefficients for sampled individuals. It is important to note that the inbreeding coefficients that we have included are exact values obtained from pedigrees. If an estimated inbreeding coefficient was used in place of the exact value, then H~ would not necessarily produce unbiased estimates in samples containing inbred individuals. H~ would also lead to a bias if relationships were misspecified. In our data example, relationships were assumed to be known, and for a data set of the size used for inferring the relationships (Rosenberg 2006) this assumption is generally sensible. However, for small data sets in which relationship inferences are uncertain, caution must be used when interpreting the bias of H~ applied to the same data from which relationships are estimated.

The estimators we have considered relate to within-population gene diversity. What if we consider the gene diversity between populations? Suppose we have samples from two populations, A and B, each containing related inbred individuals. The between-population analog of gene diversity is H^A,B=1i=1Ip^iq^i, where p^i and q^i are estimates of the frequency of allele i at a given locus in populations A and B, respectively (Nei 1987). Because the bias in within-population gene diversity estimates only arises from the quadratic p^i2 term in equation (1), E[i=1Ip^iq^i]=i=1Ipiqi (Nei 1987, p. 222), and H^A,B continues to be an unbiased estimator for between-population gene diversity in samples containing relatives.

Acknowledgments

We thank Ivana Jankovic, Yi-Ju Li, and two anonymous reviewers for helpful comments. This work was supported by NIH training grant T32 GM070449, NIH grant R01 GM081441, and grants from the Burroughs Welcome Fund and the Alfred P. Sloan Foundation.

References

  • Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998;62:1198–1211. [PubMed]
  • Boehnke M. Allele frequency estimation from data on relatives. Am J Hum Genet. 1991;48:22–25. [PubMed]
  • Bourgain C, Hoffjan S, Nicolae R, Newman D, Steiner L, Walker K, Reynolds R, Ober C, McPeek MS. Novel case–control test in a founder population identifies P-selectin as an atopy-susceptibility locus. Am J Hum Genet. 2003;73:612–626. [PubMed]
  • Broman KW. Estimation of allele frequencies with data on sibships. Genet Epidemiol. 2001;20:307–315. [PubMed]
  • Cann HM, de Toma C, Cazes L, et al. (41 co-authors) A human genome diversity cell line panel. Science. 2002;296:261–262. [PubMed]
  • Cavalli-Sforza LL. The Human Genome Diversity Project: past, present and future. Nat Rev Genet. 2005;6:333–340. [PubMed]
  • Depaulis F, Veuille M. Neutrality tests based on the distribution of haplotypes under an infinite-site model. Mol Biol Evol. 1998;15:1788–1790. [PubMed]
  • Driscoll CA, Menotti-Raymond M, Nelson G, Goldstein D, O'Brien SJ. Genomic microsatellites as evolutionary chronometers: a test in wild cats. Genome Res. 2002;12:414–423. [PubMed]
  • Gibbs JP, Martin WT. Urbanization, technology, and the division of labor: international patterns. Am Sociol Rev. 1962;27:667–677.
  • Gini CW. 1912 Variabilita e mutabilita. Studi Economico-Giuridici della R. Universita di Cagliari 3.
  • Hoelzel AR, Fleischer RC, Campagna C, Le Boeuf BJ, Alvord G. Impact of a population bottleneck on symmetry and genetic diversity in the northern elephant seal. J Evol Biol. 2002;15:567–575.
  • Jacquard A. The genetic structure of populations. New York: Springer; 1974.
  • Lange K. Mathematical and statistical methods for genetic analysis. 2nd ed. New York: Springer; 2002.
  • Li CC, Horvitz DG. Some methods of estimating the inbreeding coefficient. Am J Hum Genet. 1953;5:107–117. [PubMed]
  • McPeek MS, Wu X, Ober C. Best linear unbiased allele-frequency estimation in complex pedigrees. Biometrics. 2004;60:359–367. [PubMed]
  • Nei M. Analysis of gene diversity in subdivided populations. Proc Natl Acad Sci USA. 1973;70:3321–3323. [PubMed]
  • Nei M. Molecular evolutionary genetics. New York: Columbia University Press; 1987.
  • Nei M, Roychoudhury AK. Sampling variances of heterozygosity and genetic distance. Genetics. 1974;76:379–390. [PubMed]
  • Ohta T. Linkage disequilibrium between amino acid sites in immunoglobulin genes and other multigene families. Genet Res. 1980;36:181–197. [PubMed]
  • Ramachandran S, Deshpande O, Roseman CC, Rosenberg NA, Feldman MW, Cavalli-Sforza LL. Support from the relationship of genetic and geographic distance in human populations for a serial founder effect originating in Africa. Proc Natl Acad Sci USA. 2005;102:15942–15947. [PubMed]
  • Rosenberg NA. Standardized subsets of the HGDP-CEPH Human Genome Diversity Cell Line Panel, accounting for atypical and duplicated samples and pairs of close relatives. Ann Hum Genet. 2006;70:841–847. [PubMed]
  • Rosenberg NA, Mahajan S, Ramachandran S, Zhao C, Pritchard JK, Feldman MW. Clines, clusters, and the effect of study design on the inference of human population structure. PLoS Genet. 2005;1:660–671. [PubMed]
  • Sabatti C, Risch N. Homozygosity and linkage disequilibrium. Genetics. 2002;160:1707–1719. [PubMed]
  • Sabeti PC, Reich DE, Higgins JM, et al. (17 co-authors) Detecting recent positive selection in the human genome from haplotype structure. Nature. 2002;419:832–837. [PubMed]
  • Shete S. Uniformly minimum variance unbiased estimation of gene diversity. J Hered. 2003;94:421–424. [PubMed]
  • Simpson EH. Measurement of diversity. Nature. 1949;163:688–688.
  • Watterson GA. The homozygosity test of neutrality. Genetics. 1978;88:405–417. [PubMed]
  • Weir BS. Sampling properties of gene diversity. In: Brown AHD, Clegg MT, Kahler AL, Weir BS, editors. Plant population genetics, breeding and genetic resources. Sinauer Associates; 1989. pp. 23–42.
  • Weir BS. Genetic data analysis II. Sunderland (MA): Sinauer Associates; 1996.
  • Weir BS, Anderson AD, Hepler AB. Genetic relatedness analysis: modern data and new challenges. Nat Rev Genet. 2006;7:771–780. [PubMed]

Articles from Molecular Biology and Evolution are provided here courtesy of Oxford University Press