REFINING GENETICALLY INFERRED RELATIONSHIPS USING TREELET COVARIANCE SMOOTHING¹

GWAS panels

The human genome contains many millions of single nucleotide polymorphisms (SNPs) and other genetic variation distributed across the genome. In a GWAS it is now typical to measure a panel of at least 500,000 SNPs from each subject. SNPs typically have only two forms or alleles within a population. Whichever allele is less frequent is called the minor allele. The genotype of an individual at a SNP can then be coded as 0, 1 or 2 depending on the number of minor alleles the individual has at that SNP. Alleles at SNPs in close physical proximity are often highly correlated (i.e., in linkage disequilibrium). When multiple SNPs are in linkage disequilibrium, we say one of these SNPs “tags,” or represents, the others. Although estimates vary, well-designed panels of 500,000 SNPs do not tag all of the common SNPs in the genome and they tag very few of the SNPs with rare minor alleles [Yang et al. (2010a)]. Nevertheless, GWAS provide considerable information about familial relationships.

Estimating genetic relationships

The relatedness between a pair of individuals is defined by the frequency by which they share alleles identical by descent (IBD). Formally, two alleles are considered IBD if they descended from a common ancestor without an intermediate mutation. Within a pedigree relatives share very recent common ancestors, hence, many alleles are IBD. For a more detailed exposition of genetic relationships, see Astle and Balding (2009).

The quantity of interest in this investigation is the Additive Genetic Relationship which is defined as the expected proportion of alleles IBD for a pair of individuals. For individuals i and j we use A_ij to denote this quantity, which is more familiar when viewed as the degree of relationship, where R_ij = −log₂(A_ij). For example, for siblings, first cousins and second cousins, who are 1st, 3rd and 5th degree relatives, A is 1/2, 1/8 and 1/32, respectively. Within a noninbred pedigree A can be computed using a recursive algorithm [Thompson (1986)]. For example, if individual i has parents k and l, then A_ij = A_ji = 1/2(A_jk + A_jl).

For distantly related individuals, detailed pedigree information is not often available; however, with GWAS data one can calculate genome-average relatedness directly [Astle and Balding (2009)]. Even with complete information regarding IBD status of the chromosomes, the fraction of genetic material shared by relatives will differ slightly from the expectation calculated from the pedigree due to the stochastic nature of the meiotic process [Weir, Anderson and Hepler (2006)]. For the purpose of genetic investigations, one could argue that genome-average relatedness is a truer measure of relatedness. For example, while two distantly related individuals are expected to share a small fraction of their genetic material, if they do not inherit anything from their common ancestor, it seems appropriate to consider them unrelated.

Under many population genetic models A_ij can also be interpreted as a correlation coefficient. Let z_ik denote the scaled minor allele count for individual i at SNP k: $z_{ik}=(z_{ik}^{\ast }-2pk)/{(2pk(1-pk))}^{1/2}$, where $z_{ik}^{\ast }$ is the minor allele count and p_k is the minor allele frequency. For individuals i and j at genetic variant k, it follows from our model that

Exploiting this feature leads to a method of moments estimate of A from a panel of m genetic markers. To see this, let z_k denote a column vector of observed scaled allele counts for all individuals at the kth SNP, then let

where Z = (z₁, …, z_m). The Genome-wide Complex Trait Analysis (GCTA) software from Yang et al. (2010b) computes this estimator.

The method of moments estimator is unbiased if the population allele frequencies are known [Milligan (2003)]. In practice, the p_k’s are estimated from the sample data. A criticism of this estimator is that some off-diagonal elements are negative, which does not conform to the interpretation of A_ij as a probability. Viewed as a correlation coefficient, however, negative quantities suggest the pair of individuals share fewer alleles than expected given the allele frequencies. Alternatively, maximum likelihood estimators of A have been developed [Thompson (1975), Milligan (2003)], but these estimators are quite computationally intensive for GWAS panels. Hence, while method of moments estimators are typically less precise than maximum likelihood estimators, they are more commonly used when a large SNP panel is available.

Estimating heritability

By definition, the heritability of a quantitative trait (y) such as height is determined by the additive effect of many genes and genetic variants (g), each of small effect (i.e., the polygenic model). For individuals i = 1, …, n, suppose that the genetic effects are explained by J causal SNPs, and we can express the genetic effect as

where u_j is the additive random effect of the j th causal variant, weighted by the scaled number z_ij of minor alleles at this variant. Let g = (g₁, …, g_n)^t be the vector of random effects corresponding to the additive genetic effects for individuals i = 1, …, n. For u = (u₁, …, u_J)^t and Z_c = [z_ij], we write g = Z_cu. Define G as the variance–covariance matrix of g. Assuming $\text{Var}[\mathbf{u}]=I{\sigma }_u^2$, it follows that

where ${\sigma }_g^2=J{\sigma }_u^2$.

In the traditional model for quantitative traits a continuous phenotype y is modeled as

where e = (e₁, …, e_n)^t is the vector of residual effects, and y = (y₁, …, y_n)^t is the vector of phenotypes. In matrix notation, y = 1μ + g + e. The residuals are assumed to be independent with variance–covariance equal to $I{\sigma }_e^2$ and the random effects and residual error are assumed to be normally distributed. Consequently,

The heritability of the phenotype y is defined as

This quantity is more accurately known as the additive or narrow-sense heritability, in contrast to the broad-sense heritability, which includes nonadditive genetic effects such as gene–gene interactions. Our inferences will be confined to narrowsense heritability.

If the causal SNPs (or good tag SNPs) and the phenotype were directly measured, then one could estimate h² based on equation (5) and the implied random effects model using maximum likelihood (REML) [Searle, Casella and McCulloch (1992)]. Notationally, Z_c is an n × J matrix that picks out J columns of the full SNP panel Z. In practice, Z_c is not known. Few of the causal SNPs are known for any phenotype, and many causal SNPs will be missing from Z (i.e., not tagged by any measured SNPs).

How then is h² estimated in practice? Assuming various subsets of individuals in the sample are related with relationship matrix A (defined previously), heritability can be estimated without any knowledge of causal genetic variants that constitute g. From equation (1) and the polygenic model it follows that $\frac{Z_c{Z}_c^t}{J}\to A$ as J gets large. This inspires an alternative random effects model which has long been utilized in population genetics:

Historically, A has been derived from known pedigree structure. However, provided some subsets of the individuals in the sample are related (even distantly), one can estimate A from genetic markers using either method of moments or maximum likelihood estimation techniques. This approach has been applied frequently in quantitative genetics, especially in breeding studies [Lynch and Ritland (1999), Eding et al. (2001), Visscher et al. (2006), Hayes and Goddard (2008)]. We conjecture that by using TCS, we can improve estimates of A and obtain better estimates of heritability without knowledge of causal variants.

Alternatively, if the sample is completely unrelated, then substituting the result of equation (2) for (6) does not lead to an estimate of h² unless all of the causal SNPs have been recorded. Instead this approach estimates $h_s^2\le h^2$, the proportion of the variance in phenotype explained by the SNP panel [Yang et al. (2010a)]. In this setting, TCS will not improve estimates of $h_s^2$.

Treelet covariance smoothing (TCS)

The genetic relationship matrix A is a measure of the additive covariance structure that exists between individuals due to a common genetic background. We estimate the relationship matrix using genotyped SNPs, but this estimate is usually noisy. Hence, we propose a method for improving upon this estimate using treelets.

Treelets simultaneously return a hierarchical tree and orthonormal basis functions supported on nested clusters in the tree—both reflect the underlying structure of the data. Here we extend the original treelet framework [Lee and Nadler (2007), Lee, Nadler and Wasserman (2008)] for smoothing one-dimensional signals and functions, to a new means of smoothing and denoising variance–covariance matrices with hierarchical block structure and unstructured noise. The main idea is to first move to a different basis representation through a series of local transformations, and then impose sparsity by thresholding the transformed covariance matrix. We refer to the approach as Treelet Covariance Smoothing (TCS). The general setup is as follows. [See Appendix in Lee, Nadler and Wasserman (2008) for details on how to compute the treelet transformation. The treelet algorithm, as well as its implementation, is available in R on CRAN as the treelet library.]

Let z be a random vector in ℝ^N with variance–covariance matrix Σ. In our context, z represents the scaled minor allele counts for a set of N individuals at any SNP, and the covariance Σ = A, the additive genetic relationship matrix of the N individuals [equation (1)]. Now at each level of the treelet algorithm, we have an orthonormal multiscale basis. Let {v₁, …, v_N} denote the basis at the top of the tree [corresponding to level ℓ = N − 1 if using the notation in Lee, Nadler and Wasserman (2008)]. We write

where c_i = 〈z, v_i〉 represent the orthogonal projections onto local basis vectors on different scales. It follows that the covariance of z can be written in terms of a multi-scale matrix decomposition

where γ_i,i = Var(c_i) and γ_i,j = Cov(c_i, c_j). The first term in equation (9) describes the diagonally symmetric block structure of the variance–covariance matrix. These blocks are organized in a hierarchical tree. The second term describes a more complex structure, including off-diagonal rectangular blocks, which are also hierarchically related to each other in a multi-scale matrix decomposition.

In practice, the covariance Σ is unknown, and both the covariance matrix and the treelet basis need to be estimated from data. For relationship matrices, one can, for example, derive an estimate Σ̂ = Â from marker data using method of moments or maximum likelihood methods. Denote the treelet basis derived from Σ̂ by {v̂₁, …, v̂_N}, and write

where ${\widehat{\gamma }}_{i,i}=\widehat{\text{Var}}(c_i)$ and ${\widehat{\gamma }}_{i,j}=\widehat{\text{Cov}}(c_i,c_j)$.

Let T(Σ̂) be the covariance estimate after a treelet transformation, that is, after applying a full set of N − 1 Jacobi rotations of pairs of correlated variables. A calculation shows that

where c_i ≡ 〈z, v_i〉 and c_j ≡ 〈z, v_j 〉. This suggests² a smoothed estimate of the covariance by thresholding:

with the thresholding function

where λ is a smoothing parameter.

To summarize and in matrix short-hand notation, the smoothed genetic relationship matrix is given by

where B = (v̂₁, …, v̂_N) and T(Â), respectively, denote the treelet basis and the covariance matrix at the top of the tree, and f_λ corresponds to element-wise thresholding [equation (12)]. Note that to compute B we only need to know the Jacobi rotations at each level of the tree, more precisely, the treelet basis, B = J⁽¹⁾ · J⁽²⁾ ····· J^(N−1), where the Jacobi rotation matrix J^(ℓ) is the rotation matrix at level ℓ. The covariance estimate after a treelet transformation and before smoothing is Σ̂^ℓ ≡ T(Â) = B^t ÂB.

Choosing a smoothing parameter

The goal is to choose a threshold (λ) that reduces noise in the estimated relationships. Traditional cross-validation is not an option because we cannot predict A_ij without including persons i and j. Alternatively, we have an abundance of genetic information from which to estimate Â. We propose a SNP subsampling procedure to estimate the tuning parameter.

We begin by breaking the genome into independent training and test sets by randomly placing half the chromosomes into each set. To improve the efficiency of our estimate of A, we utilize a “blackout window” of length b to avoid sampling SNPs that are highly correlated. This b can be considered either in terms of physical location along the chromosome or the number of SNPs between any two SNPs in question. From the set of training chromosomes, select a relatively large sample of M independent SNPs to get a reliable estimate of Â. We train our algorithm by smoothing  using TCS to get Ã(λ), for all λ ∈ Λ, where Λ is a grid of reasonable threshold values.

Once we have Ã(λ), for a given λ, we subsample L SNP sets of size k from the test set of chromosomes. Here, k ≪ M and the SNPs within each of the L subsampled sets follow our defined blackout window, b. Then, for all l = 1, …, L, estimate the relationship matrix, Â_l, based on the subset of SNPs. We then compare our smoothed relationship matrix, Ã(λ), from the training chromosomes to each of the L nonsmoothed relationship matrices, Â_l, via a weighted risk function:

where w_ij is a weight associated with each element in A. Clearly, the optimal tuning parameter is λ̂ = argmin_λ∈Λ H(λ).

The reason for introducing the weighting scheme is because many subjects are nearly unrelated. Thus, we aim to upweight the loss function so that the preponderance of near-zero elements in the off-diagonal do not overwhelm the loss function. We suggest using the learned hierarchical tree to get the weights. More specifically, w_ij = |[T(Â)]_ij|, corresponding to the absolute value of the correlations between the final groupings of individuals after N − 1 rotations [equation (10)]. Also, we set w_ii = 0 because we are not interested in estimating inbreeding coefficients. It should be noted that this is a rather general weighting method. Other schemes may be more appropriate if there is a priori information suggesting the importance of particular relationships.

Simulations

To produce realistic simulations, we started with the phased genomes (haplotypes) of individuals from the HapMap 3 database³, selecting two populations with European ancestry (CEU and TSI). Utilizing the small sample of available haplotypes, our first objective was to generate a large sample of haplotypes, representative of those that might be sampled from unrelated founders of a population. The challenge was to keep intact the realistic haplotype structure for a human population, including linkage disequilibrium (LD), without generating unusual sharing between the founders. To accomplish this goal, we took the HapMap data on CEU and Tuscan samples, which were phased quite accurately into haplotypes, as the initial sample of chromosomes from which to generate founders. Now each founder haplotype was created by sampling pieces of chromosomes (or haplotypes) from the initial sample. To do so, the number of recombination spots per chromosome was determined using an overall recombination of θ = 10⁻⁶ per Mb, which is 100 times the normal rate of recombination for humans. The actual location of the recombination spots were then determined using the recombination map provided by HapMap, a procedure that successfully keeps intact the LD structure of the chromosome. From this pool of generated haplotype pairs, chromosomes were randomly assigned to each of the 39 founders in each of 100 families. These founder chromosomes were then dropped through a seven generation pedigree; see Figure 1 for the pedigree of a single family used for simulations. At each generation the chromosomes underwent recombination with an overall rate of θ = 10⁻⁸ at locations determined by HapMap’s recombination map. Within each pedigree, the genotype information of twenty individuals was collected (colored yellow). We then sampled ten individuals of varying relatedness from this group with a random sampling strategy that favored individuals of distant relatedness within the pedigree. The highest pairwise relatedness within a family is 0.125, corresponding to R = 3, and the lowest is < 0.001. Individuals from different families are unrelated. Each simulation produced a total of 1000 individuals made up of 100 ten-member families of varying levels of relatedness. Finally, the entire process was repeated fifty times.

Because we know the pedigree structure, we can compare the unsmoothed estimate  to à found via TCS. Here, we use the GCTA software [Yang et al. (2010b)] to estimate  using 100,000 randomly chosen SNPs with MAF > 0.05. The optimal level of smoothing (λ̂) is chosen via the subsampling scheme described previously using M = 5000, b = 10, k = 50, L = 50 and repeating everything ten times. Here, b is in terms of number of SNPs. We choose λ̂ by examining a plot of H(λ) across a grid of λ values. The optimal smoothing parameter is the one that minimizes the risk function, H. For one such simulation sample we can see from Figure 2 that λ̂ ≈ 0.051.

The question then becomes, does TCS improve estimates of relatedness? Figures 3 and 4 display boxplots comparing the root mean square error (RMSE) of  to à at varying levels of known pairwise relationship values. For a full comparison, we have included two smoothing methods: TCS, as previously described, as well as “simple thresholding,” wherein the elements of  are directly thresholded. [The latter approach is a degenerate case of TCS models, at level ℓ = 0, for which the basis is the Dirac basis, i.e., v_i = v̂_i = δ_i for i = 1,…, N in equations (8)–(13).] Moving from left to right in the figures, the true degree of relatedness increases from R = 4, …, 11, to no relatedness. Over the entire matrix of estimates, the RMSE is 0.0055, 0.0015 and 0.0019 for the unsmoothed (Â), TCS (Ã_t) and simple thresholding (Ã_s) methods respectively, demonstrating an overall advantage of TCS. As with many shrinkage methods, TCS introduces a slight bias that is reflected in a higher RMSE for closely related individuals. Consequently, TCS has a larger RMSE than the unsmoothed estimate for smaller values of R. Where TCS gains a notable advantage over the unsmoothed estimate is in differentiating between more distantly related individuals and noise. From Figures 3 and 5 we can see that simple thresholding incurs a substantially larger RMSE for closer relationships because it thresholds too aggressively. For R = 4, 70% of the pairs are zeroed out, and for R > 4 virtually all pairs are zeroed out. Naturally, this method has the smallest RMSE for the sample of unrelated pairs because thresholding zeros out all of these entries. Notably, TCS does almost as well in this setting. A direct comparison of RMSE does not fully reflect the true loss incurred in practice. In most genetic studies close relatives are often recorded in pedigrees and, hence, estimates are not required. Alternatively, considering distant relatives to be unrelated leads to a substantial loss for estimating heritability and most other genetic applications.

Heritability in health ABC study

Body Mass Index (BMI) is one of several traits measured as part of the study entitled “Whole Genome Association Study of Visceral Adiposity” as part of the Health Aging and Body Composition (Health ABC) Study. These data are archived in the Database for Genotypes and Phenotypes (dbGaP)⁴. We restrict our attention to those 1644 individuals with self-reported European ancestry. To control for confounding, prior to analysis, we adjust BMI scores by regressing out age, gender and collection site. Our objective is to estimate heritability of BMI from this population sample. Published heritability estimates range from as low as 0.05 to as high as 0.90 [Allison et al. (1996)]; however, based on estimates derived from known pedigrees, the heritability of BMI is estimated to be approximately 50–75% [Kangas-Kontio et al. (2010), Zabaneh et al. (2009)].

From the full sample of SNPs (Illumina 1M platform) we remove those with missingness greater than 0.1% and MAF < 0.01. From these we select a subpanel of 90,000 SNPs, chosen to be nearly evenly spaced. Based on these SNPs, we calculated the relationship matrix Â, and find that the individuals are predominately unrelated. The most highly related pair appear to be third degree relatives, such as first cousins. And more than half of the pairs appear to be more distantly related than 10th degree relatives.

To estimate the heritability in this setting, we input the smoothed relationship matrix in equation (11) into the REML algorithm. The required smoothing parameter λ is selected in two ways: (i) minimizing the loss function in equation (14) via the subsampling approach; and (ii) using a profile likelihood approach. With both techniques, we get estimates of the heritability that are very close to what is found in the literature.

For a range of smoothing parameters, 0 ≤ λ ≤ 0.40, we calculate the smoothed relationship matrix, Ã_λ, and plug this value into the REML model to obtain a profile likelihood (Figure 6). Also plotted in this figure is ${\widehat{h}}_{\lambda }^2$, the heritability that maximizes REML as a function of λ (or minimizes—2 times the log-likelihood). Without smoothing (λ = 0), which is not shown in the plot, ĥ² = 0.23. Smoothing the relationship matrix results in an increasing estimate of the heritability which stabilizes at about 70%. Further smoothing beyond the range displayed leads to a numerically unstable optimization problem and diminished likelihood. Using the profile likelihood approach, λ is chosen to be the point at which REML is maximized. This method results in an estimate of λ̂ = 0.20 corresponding to ĥ² = 0.71. Smoothing using our SNP subsampling scheme results in λ̂ = 0.18 and ĥ² = 0.72.

For comparison, we have repeated the above experiments with an orthogonal basis computed by principal component analysis (PCA) in lieu of a treelet basis. Such an approach does not improve the estimates of family relationships or heritability. When noise is present, PCA is unable to uncover the underlying sparse structure of the relationship matrix. In fact, the results with PCA are identical to those without smoothing (with the profile likelihood peaking when the tuning parameter is set to 0).

Another trait that was measured in this study is the Abdomen Visceral Fat Density (AVFD). As was the case with BMI, we restricted our attention to individuals of European descent and regressed out age, sex and collection site. According to the literature, the heritability of AVFD should be between 45–70% [Katzmarzyk, Perusse and Bouchard (1999)]. According to Figure 6, one can see that using the smoothing parameter based on our subsampling scheme (λ̂ = 0.18) we get ĥ² = 0.29. On the other hand, exploiting the profile likelihood plot results in λ̂ = 0.09 and ĥ² = 0.36. When no smoothing was used (not shown in figure), λ̂ = 0.11. Thus, both methods for choosing the smoothing parameter used in TCS resulted in estimates of the heritability that are closer to what is established in the literature than without smoothing.

It is notable that ĥ² for both traits increased toward the established estimate of heritability regardless of how we estimate the optimal smoothing parameter, because only a small fraction of the genome was sampled by the SNP panel. Thus, our results underscore the fact that the quantitative trait model given in equation (5) does not require measurement of the causal SNPs that constitute equation (3). What is required is a good estimate of A based on relatives.

Our analysis of BMI and AVFD illustrates the difference between estimates of heritability in the traditional sense and estimates of $h_S^2$, the heritability attributable to the SNPs in the panel. From equations (6) and (7) it is clear that heritability derived from the classic quantitative traits model can distinguish between variance explained by relatives and variance explained by causal SNPs only if either (i) all causal SNPs are excluded, or (ii) all relatives are excluded. Because a large number of undiscovered SNPs scattered across the genome are likely to be causal, and large samples invariably contain distantly related individuals, some ambiguity will always be present.

Clearly, the 90,000 SNPs in our panel do not explain a substantial fraction of the variation in BMI and yet we obtain an accurate estimate of heritability using TCS. The increase in estimated heritability of BMI from 23% to 72% suggests that smoothing improves the estimate of A and that a substantial fraction of the correlation in BMI in our sample is due to genetic relatedness. In a similar study with a larger population sample Yang et al. (2011) estimated $h_S^2$ of BMI at 17% when using the full SNP panel, but excluding all detectable relatives. Assuming relatives were successfully removed, they conclude that approximately 17% of the variability in BMI is explained by common variants included or tagged by the SNP panel.