Introduction

Most of the complex genetic diseases have risk factors that are quantitative in nature. For instance, cholesterol level is a risk factor for myocardial infarction, Body mass index is a risk factor for type 2 diabetes. These quantitative traits (QTs) often have strong genetic determinants. It is therefore, of considerable interest to map the genes underlying QTs [1]. For most QTs relevant for human health and disease a large proportion–ranging from 30 to 80%–of variation is explained by genetic factors. Multiple genes are expected to contribute to this variation. The proportion of variation explained by a single gene is expected, however, to be small (less than 5%). For example, one of most known quantitative trait loci (QTLs), APOE is strongly and consistently associated with increased total cholesterol levels. Yet it explains only about 2–5% of the variation of this trait [2], [3]. When effects of specific common variants are expected to be small, association analysis provides a powerful approach to identify the gene compared to linkage analysis [4]. Genome-wide association analysis is a powerful tool to disentangle the complexity of quantitative traits, even in family based studies.

For pedigree-based association analysis several methods and software packages have been developed that utilize information about transmission of alleles, such as the orthogonal test for within-family variation [5] and family-based association test [6], [7]. As these methods analyze within-family variation, they are robust to population stratification. However, these methods ignore a large proportion of information provided by the between-family variation leaving room for improvement.

A conventional polygenic model of inheritance, which is a statistical genetics' “gold standard”, is a mixed modelwhere μ is the overall mean, G is the vector of random polygenic effect, and e is the vector of random residuals. This model may be extended to study association by including a kg term, where k is the marker genotype effect, and g is the vector of genotypes. Such a model is known as the measured genotype (MG) model [8]. The MG approach, implemented using (restricted) maximum likelihood, is a powerful tool for the analysis of QTs when ethnic stratification can be ignored [9], [10] and pedigrees are small or when there are few dozens or hundreds of candidate polymorphisms to be tested. This approach, however, is not efficient in terms of computation time. This hampers the application of MG in genome-wide association analysis.

We proposed a fast two-step approximation to MG, a Genome Wide Rapid Association using Mixed Model and Regression (GRAMMAR) [11] . In the first step the individual environmental residuals are estimated, using additive polygenic model. Then the test of association is performed using these familial correlation-free residuals with a rapid least squares or score method. Though the two-step method is indeed computationally efficient and outperforms family based approaches like FBAT and QTDT in terms of power and speed, it loses power compared to MG [11]. The test becomes increasingly conservative and less powerful with the increase in the number of large full-sib families and increased heritability of the trait. Interestingly, empirical power of GRAMMAR is very close to that of MG.

Both classical MG and GRAMMAR approaches rely heavily on the availability, completeness, and correctness of genealogical information. When these assumptions are violated, the most likely outcome is inflation in type 1 error. Practically, genealogical information may often be available only for a limited number of generations, it may be inaccurate, and it may become increasingly inaccurate back in time. This may be taken as an argument for application of underpowered TDT-like methods to avoid false positive or negative results.

We and others reason that genomic data can be used to correct for the (only partly observed) true genealogy. With the new array technologies, large numbers of markers can be typed over the genome. These provide information on ‘genomic background’, which can be used to infer population (sub)structure and relations between study participants, which is classically done using the genealogy. In a recent study of type 2 diabetes, genomic control (GC) was applied to control for relatedness among cases and controls from Icelandic population [12], However, the type 1 and 2 errors of GC were not yet systematically investigated in the context of pedigree data analysis. In this study we aim to exploit the ideas of genomic background to extend our work on family based association [11] and determine how powerful and efficient the method for genome-wide association analysis of QTs in samples of related individuals is.

Methods

Only minor proportion of markers in a genome-wide association study is expected to be truly associated with an analysis trait, and a vast majority of the markers may be thought of as realizations under the null hypothesis and can be used to characterize the null distribution of the test statistics. This idea follows that of Genomic Control (GC) method [13], which was introduced in the context of association analysis in population-based studies, where population stratification and cryptic relatedness may be present.

Following Devlin and Roeder [13], [14] we suggest estimation of the test statistic inflation factor λ by regressing the trait on N loci. From each regression analysis, estimate T_i² = β_i²/Var(β_i), where β_i is the effect of the i-th SNP (i from 1 to N) and Var(β_i) is the variance of the estimate. Inflation factor is estimated aswhere 0.456 is the median of the χ₁² distribution with a non-central variance equal to φ. The number of loci used, N, in a genome-wide association study is typically reflecting all loci investigated or the ones generating lowest 95% of test statistics [15], [16]. It should be noted that the value of λ in conventional GC is constrained to be greater than one as values less than one have no theoretical meaning.

We propose use of GC to correct for conservativeness of the GRAMMAR approach outlined earlier. This method which we call GRAMMAR-GC involves three steps: (a) trait heritability is estimated by using trait and pedigree data using the following mixed modelwhere β_j is the effect of j^th covariate , c_ij is the value of j^th covariate and μ, G, and e are defined earlier. And environmental residuals are estimated as(b) the markers are tested for association with these residuals using simple linear regressionwhere k and g are defined earlier

(c) GC procedure is applied to correct the test statistic. The deflation factor ζ is estimated by regressing residuals from step (a) on each of the k null loci and from each regression analysis is estimated, where is the effect of the k^th SNP. The deflation factor ζ is estimated as . Then T²/ is compared with χ²₍₁₎ to determine whether the locus is significantly associated with the QT [14].

Steps (a) and (b) comprise the original GRAMMAR approach, leading to a conservative test. In step (c), GC is used to estimate the deflation factor ζ. This deflation factor is estimated in exactly the same way Bacanu et al. [14] estimate inflation factor λ for quantitative traits. The only difference is that ζ<1 in contrast to λ in conventional GC methods which is constrained to be >1. This difference comes from the fact that residuals from step (a) are regressed on null loci to obtain the estimate of ζ instead of original trait data as in conventional GC.

The original GRAMMAR relies on the availability of a precise and complete pedigree structure for heritablity estimation in step (a). This can, however, be done by using kinship coefficients estimated from genomic data where the genomic estimate of kinship for a pair of individuals i and j is obtained using the formula [17]where g_ik is the genotype of the i-th person at the k-th SNP of the (coded as 0, 1/2 and 1, for rare allele homozygote, heterozygote and common homozygote, respectively), p_k is the frequency of the major allele, and n is the number of SNPs used for kinship estimation.

Heritability is then estimated by maximizing the log-likelihood of the datawhere y is the vector of trait values, μ is the mean and Σ = Φσ_G²+Iσ_e² is the variance-covariance matrix. Here, Φ is the relationship matrix whose elements are defined as 2 f_ij, σ_G² is the additive genetic variance due to polygenes, I is the identity matrix and σ_e² is the residual variance. Trait environmental residuals were obtained as

To avoid confusion we refer to the new method as Pedigree GRAMMAR-GC (PedGR-GC) when in step (a) heritability is estimated from the genealogy, and Genomic GRAMMAR-GC (GenGR-GC) when the heritability is estimated from the genomic data. For PedGR-GC, environmental residuals were estimated using ASReml [18]. All other computations were performed using freely available R software (http://www.r-project.org); computations associated with GenGR-GC were facilitated by GenABEL package for R [19], implementing procedures to compute genomic kinship, maximize polygenic models and compute the residuals.

Results

Table 1 shows the 95^th percentile of the test statistic and type 1 error (proportion of simulations that resulted in a χ²≥3.84 and χ²≥6.64, corresponding to tabulated α = 0.05 and 0.01 respectively) for GC and PedGR-GC and GRAMMAR. For GC and PedGR-GC type 1 error is close to the nominal α while GRAMMAR is conservative and this conservativeness increases with the increase in the number of sibships and the heritability of the trait.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. 95^th percentiles of the distribution and type 1 error for PedGR-GC, GC and GRAMMAR.

https://doi.org/10.1371/journal.pone.0001274.t001

Supplementary Table S1 and Figure 1 illustrate the power of the proposed (GC and PedGR-GC) and previous methods (MG and GRAMMAR). In the Figure 1, power to reach χ²≥6.64 (α = 0.01) is depicted with circles. From the available evaluation points we also estimated the slope of linear regression of non-centrality parameter on proportion of variance explained and used this slope to derive power curves.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Power of MG (red line), GRAMMAR (green line), PedGR-GC (blue dashed line), and GC (pink dashed line) to detect association under different heritability models and pedigree structures.

The three rows show the power under different heritability models (from 30% to 80%) and the three columns show power achieved in different pedigrees namely nuclear pedigrees (NP), Erasmus Rucphen Family (ERF), and idealized pig population (IPP). The y-axis of each panel shows power while the x-axis shows the proportion of variance explained by the QTL under study. The red (for MG), green (for GRAMMAR), blue (for PedGR-GC), and pink (for GC) circles show the empirical power estimates. The power estimates are based on α = 0.01. The empirical power estimates are based on 1000 simulations for NP, and IPP, and 100 simulations for ERF.

https://doi.org/10.1371/journal.pone.0001274.g001

Figure 1 shows that the power of PedGR-GC (blue dashed line) is very close to the power of the ‘gold standard’ measured genotype approach (red line) for all scenarios. These two methods appear to be the most powerful of all methods for all pedigree structures and genetic models evaluated. The power of GC (pink dashed line) is close to that of MG and PedGR-GC when the heritability is low but its' power declines when the heritability of the trait increases or when pedigrees with large number of full-sib families (IPP) are investigated. GRAMMAR (green line) performs similar to GC in a sample of nuclear families and in the ERF sample, but is more powerful when IPP pedigree is investigated.

To study the potential of described methods on genome-wide scale we used 695 people who are a part of the ERF pedigree and who were genotyped using Illumina 6K SNP linkage panel. Based on pedigree records, the 695 people formed 471 pairs of first-, 311 pairs of second-, 681 pairs of third- and 1,105 pairs of forth-degree relatives and 223,578 pairs of more distant relationship.

We generated 500 replicas for each of the models assuming total trait heritability of 30, 50 and 80%. In each replica, we selected five hundred random SNPs each explaining equal proportion of variance amounting to the total heritability minus 4%, and one random SNP explaining 4% of the phenotypic variance. The analysis trait was obtained as a sum of the SNP effects and a random number from the Normal distribution with mean zero and variance 0.7, 0.5 or 0.2 for total trait heritabilities of 30, 50 and 80%, respectively.

Type one error was estimated as the proportion of non-associated SNPs (>2.5 million tests in total) showing association P-value of 5% or less. The statistical power was estimated as the proportion of replicas in which the SNP explaining 4% of variance passed genome-wide significant threshold (Bonferroni-corrected P-value 0.05/5524 = 9×10⁻⁶).

For analysis, we used GC, PedGR-GC and GenGR-GC. For GenGR-GC, the kinship matrix used was estimated from genomic information on 5524 autosomal SNPs typed in 695 members of the ERF study.

All methods showed a genome-wide type 1 error which was very close, but lower than the pre-specified threshold of 5% (Figure 2A). The methods tended to be more conservative at higher heritabilities. These results are consistent with the observations of others, that GC is conservative, and can be explained by the fact that all SNPs, some of which were associated with the trait, are used to estimate the null distribution of the test statistic.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Type 1 error (A) and power (B) to achieve 5% genome-wide significance at the truly associated SNP in a study of 695 ERF people genotyped for 5524 autosomal SNPs.

https://doi.org/10.1371/journal.pone.0001274.g002

Discussion

In this work we aimed to develop fast and powerful methods for genome-wide association analysis in samples of related individuals by exploiting the ideas of genomic background for correction of the distribution of the test statistics and for inferring the relation between study subjects. We show that methods, which exploit only genomic background, such as Genomic Control (GC) and GRAMMAR-GC using genomic kinship, are powerful and genome-wide feasible methods. Moreover, genomic GRAMMAR-GC, which infers genetic relations from genomic data, may be superior to the methods that use pedigree information in an exact manner.

Our simulations show that GC is a valid method to study data coming from samples of related individuals. GC can be a powerful tool for the analysis of pedigree based quantitative trait loci. It outperforms traditional family based approaches like QTDT and FBAT (cf. Table 1 & Table 2 of Aulchenko et. al.[11]). The power of GC is close to that of the ‘gold standard’ measured genotype approach when trait heritability is low or moderate and human pedigrees are studied. However, the power of GC drops notably with high trait heritability and when pedigrees involve very large sibships, which are likely to be observed in animal pedigrees.

The results that GC was less powerful, than GRAMMAR-GC based methods, and tended to lose power at higher heritabilites (Figure 2B), are not surprising and are consistent with our previous findings [11]. Interestingly enough in a study of real Genome-wide data in ERF pedigree, GenGR-GC was consistently more powerful than pedigree-based GRAMMAR-GC (PedGR-GC), and the power advantage became more pronounced at higher heritabilities.

We proposed an extension to the GRAMMAR method [11], which increases its' power, maintains a nominal type 1 error and also does not require the precise knowledge of pedigree structure. Our method (GRAMMAR-GC) involves three steps which include removing the correlation from the data using relationship matrix estimated from either the pedigree or the genome (derivation of environmental residuals), using the correlation-free residuals from step 1 as the trait to perform association analysis, and then applying GC to correct the test statistic. We show through simulations that our proposed method performs very similar to the Measured Genotype (MG) approach with respect to type 1 error and power yet it is fast and feasible for genome-wide association analysis. By analysis of real genome-wide SNP data we showed that when the genomic data is used to estimate the relationship matrix (GenGR-GC) instead of the estimate obtained from genealogy, the power might be even improved.

One of our conclusions is that in genome-wide association studies of related individuals genomic background based methods such as genomic GRAMMAR-GC should be preferred over the ones exploiting known pedigree structure, such as pedigree GRAMMAR-GC or MG approach. There are two reasons why we believe that GenGR-GC should be preferred over its' pedigree analog.

First, errors in genealogy such as mis-specification of relations can lead to either false positives or negatives. Secondly, relationship coefficient computed from a pedigree is an expectation of the proportion of genome shared identical by descent (IBD) under the infinitesimal model, assuming infinite number of unlinked loci. The true proportion of genome shared, however, may deviate from this expectation [17]. For example, for remote relatives there is a fair chance of not sharing any genomic loci IBD. We may speculate that kinship estimated based on marker data can reflect the true unobserved genomic sharing better then the expectation computed from even a correct pedigree. If this is true, under the polygenic model we should expect that methods based on kinship estimated from marker data will be more powerful than the methods estimating kinship from the pedigree. We however leave a more detailed investigation of effects of pedigree error and precision of genomic kinship estimates on type 1 error and power to future works.

Another advantage of the GRAMMAR-GC is that the environmental residuals used for the analysis are free from familial correlations. Therefore the structure of the data becomes exchangeable. This means that permutation techniques may be applied to derive empirical measures of significance. In the analysis of data where adjacent SNPs are correlated due to linkage disequilibrium, thresholds set via permutation will account for these correlations and result in less stringent thresholds than those set by Bonferroni correction. This will become more and more important in the future, when denser marker maps with millions of SNPs will be applied to do association studies. Another attractive feature is that a range of new methods developed for classical “unrelated individuals” design can be applied to polygenic residuals obtained at the first stage of the analysis. In recent years, much progress was made in development of powerful methods and software which are robust to possible allelic heterogeneity through utilization of haplotype clustering and population genetic coalescence modeling [21], [22].

Finally, GRAMMAR-GC is easily extendable: for example, it is easy to include covariates, interactions with sex and environment, gene-gene interactions and parent of origin effects.

To conclude, GRAMMAR-GC is a fast powerful approach for genome-wide association analysis of quantitative traits in samples of related individuals, which does not require precise knowledge of pedigree structure. This method is implemented as part of the GenABEL package, available at http://cran.r-project.org/.

Supporting Information

Acknowledgments

We would like to thank Prof. Tatiana Axenovich, Prof. Haley, Dr de Koning, Dr. Catherine Bourgain, Prof. Ober, and Dr. William Astle for valuable discussions.