Sequencing technologies are becoming cheap enough to apply to large numbers of study participants and promise to provide new insights into human phenotypes by bringing to light rare and previously unknown genetic variants. We develop a new framework for the analysis of sequence data that incorporates all of the major features of previously proposed approaches, including those focused on allele counts and allele burden, but is both more general and more powerful. We harness population genetic theory to provide prior information on effect sizes and to create a pooling strategy for information from rare variants. Our method, EMMPAT (Evolutionary Mixed Model for Pooled Association Testing), generates a single test per gene (substantially reducing multiple testing concerns), facilitates graphical summaries, and improves the interpretation of results by allowing calculation of attributable variance. Simulations show that, relative to previously used approaches, our method increases the power to detect genes that affect phenotype when natural selection has kept alleles with large effect sizes rare. We demonstrate our approach on a population-based re-sequencing study of association between serum triglycerides and variation in ANGPTL4.
Several authors have reviewed the potential contribution of low frequency alleles to variation in phenotypes [2]–[7]. Absent a change in the properties of new mutations during recent history, which we find implausible, systematic differences between SNPs of varying frequencies must be mediated by natural selection. Since the early 20th century, much work has explicated the evolutionary dynamics of quantitative traits, reviewed by Barton and Johnson [22], [23]. Below we will posit a model of pleiotropic selection whereby the trait under study or a trait with a correlated genetic basis is under purifying selection. More detailed connection and contrast to the existing work on the genetic basis of quantitative traits is found in Text S1. In Figure 1, we illustrate direct and apparent selection scenarios which give rise to a correlation between fitness effects and phenotype effects. In Figure 1A, the phenotype itself is under selective pressure; for example, disease leading to propensity to childhood mortality. Figure 1B shows apparent selection by pleiotropy; variants which disrupt an unconstrained role of a gene also tend to disrupt another role which is under selection; for example, variation which increases Alzheimer' s Disease risk after reproductive age may relate to other brain function which is relevant for individuals still reproducing. Hartl and Clark [24] carefully constructs and interprets the concept of fitness-effects in classical population genetics. Briefly, in an idealized population, the relative reproductive advantage of an individual is the product of the fitness effects of each variant that person carries, an additive approximation with no dominance or epistasis. We parameterize the problem in terms of the log of multiplicative fitness effects. That is, the fitness of the person is given by where the fitness effect of the variant is denoted and is the unphased genotype at that locus. The fitness effect of a new mutation determines several of its properties, such as average sojourn time before either going extinct or fixing at 100% prevalence and average frequency when sampled at a point in time [24]. Rather than assume that all variants in the region have the same, we assume that the of new mutations are sampled from a distribution of fitness effects (DFE). Just as a fixed would determine properties of the sampled genotype data for a SNP, a DFE along with mutation, recombination, and demographic parameters induces a distribution on the observed frequency spectrum and polymorphism - divergence ratios in sampled data. Several authors have attempted to fit a parameterized DFE from genomic data [25]–[34]. Boyko et al [33] found that a combination of a point mass at neutrality (not under selection) combined with a gamma distribution for deleterious differences from neutrality to be a good fit for the DFE of non-synonymous mutations. With these facts in mind, in what follows we will use fitness effects to operationalize the construct of functional status for each SNP. Whereas Johnson and Barton [23] worked directly with the joint distribution of fitness and phenotype effects, we will use an existing DFE estimate [33] as a marginal distribution for fitness effects and construct the conditional distribution of phenotype effects. Since we do not know the true fitness effects of SNPs, we will estimate them with observed SNP frequency, which is statistically ancillary to phenotype-SNP correlation, using a simulation methodology described below.
Assume the context of a simple random cross-sectional sample of individuals (indexed by) studying a quantitative trait measured once per individual. Assume that these individuals also possess vectors of covariates and genotypes at each locus inside a sequenced candidate gene or region. The genotypes are coded such that “0” represents homozygous possession of the ancestral allele, “1” heterozygosity, and “2” homozygous possession of the derived allele at the locus. That is, represents the fourth sampled person possessing two derived alleles at the third locus in the sequenced region. We can write a regression model for person' s phenotype in terms of deviation from an average level predicted by covariate effects and additive genotype effects, (1) Using standard least-squares regression to estimate such a model will pose several problems. First, because there will be many rare variants, will contain many poorly estimated coefficients. The large number of rare variants will give model (1) a large number of degrees of freedom, decreasing its power to detect association with the candidate gene. Some of the variation uncovered may be perfectly correlated in the sample, meaning that those coefficients are not separately estimable in least-squares regression. Additionally, as the amount of the genome sequenced becomes large, there will be more variants than participants, making the entire model unidentified. To overcome these problems, we need to make more assumptions and model the coefficients. We adopt a model where we view the effects of SNPs in the study as a sample from a wider population of SNP effects, and characterize that entire population using only three parameters. To fix ideas, assume for now that we knew the fitness effect of each SNP. If fitness was perfectly correlated to effect on phenotype, we would use that as a summary for all alleles, , where the parameter relates the scales of the two measures. As the fitness effect is not perfectly correlated to effect on phenotype, we add a mean and an error term acknowledging those limitations to obtain (2) In applied problems, is not known a-priori, so we will construct a prediction based on the observed frequency. We denote for that estimate, and for its prediction error we write. We plug those estimates in to (2) to obtain (3) and combine the two uncorrelated error terms to yield (4) where (5) (6) The first term in (4), , allows derived alleles to on average increase or decrease the phenotype. The second term is an unscaled correlation between phenotype effects and expected fitness effects. The error term is the deviation in SNP j' s effect on phenotype from the average of SNPs with the same observed frequency. The variance of in (6) therefore has two components, first corresponds to prediction error of, and second is the variance of phenotype effects for SNPs at the same level of true fitness burden. The function allows that as average burden changes the variability might also change. Although one could imagine “bad” alleles being more variable in their effects than relatively neutral alleles, implying non-constant, we propose constant as a reasonable modeling start. This will still allow for the variance of effect sizes to change with observed frequency because of non-uniformity of with frequency. Equation (4) asserts that phenotype-effect and fitness-effect are linearly related; that seems correct for the scenario in Figure 1A and a good starting place for the other possibilities. In future work we will be able to empirically examine this assumption by graphical diagnostics and comparing fits using other functional forms. Further discussion of nonlinear relationships is found in Text S1, and we will demonstrate the impact of an incorrect assumption of linearity in our simulation studies. Our model is quite general in that existing methods correspond to submodels of (4). An allele count method tests the model with only allowed to vary; rare alleles below an arbitrary threshold are summarized by an average effect which does not change with frequency, so and all are set to zero, and alleles above that threshold are regarded as free parameters. Similarly a weighted-burden method corresponds to the model with a particular implementation of, such as in Madsen et al [21] where, and forces all in the rare alleles to be zero. Our model will not involve an arbitrary threshold for “rare alleles” and will adaptively pool variant effects in a flexible way. As shown in the results, this will create substantial power gains in a variety of settings. When and in (4) are zero, our model reduces to a standard random-effect model identical to that of Kwee et al [19] with all variants given the same weight. That is, regardless of frequency all SNPs have the same likelihood of having large effect sizes, and regardless of frequency SNP effects have zero mean. As a result, our method will be robust to the case that fitness and phenotype effects are unrelated by estimating and retaining the flexibility of the method of Kwee et al. The major difference between the above and our method is the use of population genetics to suggest the structure of the variance of SNP effects, including a fallback should fitness and phenotype effect not be related. Kwee' s method is developed in the context of tag SNPs and suggests an arbitrary variance of SNP effects given as either a constant, , or any prior-information based form. A related method is that of Hoggart et al [18]. Their approach corresponds to and set to zero (they assume a mean-zero distribution) and a different set of restrictions on the distribution of. Their assumptions about the distribution of were chosen to yield estimates with most variants having zero effect, a feature called model selection which eliminates small effects and correlated variables. In contrast, our model will tend to reign in large effect sizes and split effect size between variants in high linkage disequilibrium, but does not eliminate SNPs from the fit. We prefer our choice for resequencing for several reasons. First, there may well be many effects of small size which are cumulatively important, and we want to retain those small effects in the model. Second, we want an estimate of the effect size of each variant for graphical and diagnostic purposes. Third, we accomplish a similar goal of reducing the model size by rejecting the null on a small number of genes. That is, we want to identify a small number of disease relevant genes with our efficient test; doing so will exclude most SNPs without further model selection procedures. Fourth, by smoothly grouping rare SNPs and summarizing them with only a few parameters, we already greatly reduce the multiple testing burden. The specification of equations (1), (4), and (6) yields a natural interpretation to the fitted model. After estimating the population parameters of phenotype effects, we will be able to jointly estimate individual SNP effects and their impact on the phenotype of each person in our sample. By calculating, we obtain the expected difference between participant i' s phenotype and what we would expect were there no effects of this gene. As a result we can empirically estimate the overall phenotypic variability due to observed genetic variants, over study participants. We can similarly estimate the variability dues to rare alleles by including only rare SNPs in the above calculation. The overall effect is an average change in phenotype per derived allele, perhaps due to inadequate purifying selection. In the variance expression (6), is the variability of allelic effects for a given level of true fitness. As will be shown below in Figure 2, when using the genome-wide distribution of fitness effects for non-synonymous SNPs, common variation is nearly neutral so can also be thought of as the variability of effects of common alleles. represents the correlation between fitness burden and phenotypic burden. This parameter' s interpretation relies on accurately estimating the scale of fitness effects and has awkward units, but we can avoid this difficulty by noting that (4) can be decomposed into a fitness related portion and a fitness unrelated portion which are independent (7) By calculating we can ascribe a proportion of total variation in phenotype to selection-phenotype correlation without worrying about having gotten the scale of correct. Calculations for separating these variance components are found in Text S1. We can use the same technique to compare classes of SNPs, for example non-coding vs missense, by jointly fitting separate and comparing the attributable variance for each class of SNPs. We will illustrate this idea in our real data example. This decomposition also shows why it is not crucial for our estimates of fitness to be perfect. The model can fall back by setting to zero and use only to recover a working model which does not pool information across rare alleles. Doing so will mean that the opportunity to gain information by recognizing structure in phenotype effects will not be realized, but the remaining estimation method is still valid. An important consideration is how to interpret the results when multiple ethnic groups are analyzed simultaneously. Because some genetic variation is fixed between ethnic groups in the sample, the average effect of single-population variation will be absorbed into the fitted mean for that group. As a result, the interpretation for “total explained variation” is actually “total explained within-ethnic-group variation; ” genetic variation may explain some of the phenotypic difference between groups, but we do not include it in our estimate because of confounding between environmental exposures and ethnic background. Another point requiring clarification is the assumption that genotype effects are independent. In the context of GWAS, nearby SNPs often are thought to have correlated effects because they mutually tag a functional variant. Additionally, estimates of SNP effects will be correlated due to LD making their true separate effects difficult or impossible to identify. However, in the underlying data generating mechanism true genotype effects are independent. Because sequencing identifies all the variation within the region and eliminates much of the correlation due to untyped alleles, we believe that the independence assumption is a useful approximation in this case. Non-independence of the true effects could be accommodated by imposing a covariance structure on SNP effects, for example using their spatial distance in the genome or folded protein. Alternatively, the phylogenetic approach of TreeLD [35] estimates the degree of probable overlap of untyped SNPs. Model (4) relies on a prediction of the fitness effect of each variant as well as an estimate of the error of that prediction. We use the following procedure to calculate such estimates. To reduce computational requirements, steps 2 and 3 above can be replaced by simulating a smaller number of large populations and calculating the expected mean and variance of fitness using simple random sampling. Figure 2 depicts the relationship of and to frequency when using a genome-wide fitted DFE [33]. Because much of the variation discovered in our multi-ethnic example dataset is confined to one ethnicity, we use the ethnicity-specific frequency and pseudo-data. Because of admixture in our sample, we use the highest observed frequency (the most skeptical about its being rare) to assign an ethnicity of origin to SNPs appearing in multiple groups. An advantage of this method is that because it refers to a feature of genetic history rather than a phenotype, it need only be done once for any trait under study on the same cohort. While the fitness - phenotype relationship will be different for all traits, that is modeled by the fitted parameter rather than modification of. If the impact of LD structure on the prediction does not vary too much between genes, the calculation can be recycled for multiple genes under study. In some experiments, we found the impact of LD to be minimal (data not shown). Discussion of taking the DFE as known versus estimating or using some other flexible function of frequency it is included in Text S1. Discussion of the quality of the existing DFE estimates are also included in Text S1. We have used the observed frequency to estimate the fitness effect, but there are many other potential predictors of functional status. Discussion of including them in our model is found in Text S1.
We propose a novel method, EMMPAT, for association between sequenced genes and phenotype which utilizes population genetic theory to pool information among rare variants. Our method generalizes allele-count and allele-burden techniques, and presents several advantages. Of greatest importance to the practicing scientist will be increased power and interpretability. As shown above, our method allows us to leverage allele frequency as auxiliary data related to SNP effects and to substantially increase power to detect association in many scenarios. The availability of a well motivated pooling strategy allows an omnibus test which incorporates common and rare variation simultaneously. Our approach provides clear interpretations for the fitted model, such as the attributable variance in phenotype due to all polymorphisms observed in a gene, particular types of SNPs, or only the rare variation. Furthermore it facilitates tests of meaningful parameters (such as mean derived allele burden) and group differences (such as non-synonymous versus non-coding). The regression toolbox allows model checking and exploration, such as in Figure 3 which presents the data in an informative format. Additional model checking proceeds as usual in linear mixed models, and posterior predictive checks are similarly possible. A relevant question is how important our method will be for diseases which have not been strongly selected against. There are three answers to consider. First, when selection and disease effect are completely independent, common SNPs will tend to have just as large effect sizes as rare SNPs and explain much of the heritable variation in phenotype [2], [3]. We believe that most investigators conducting resequencing studies assume rare variation to have larger effect sizes, since that is the best-justified scenario for the expense of sequencing. Second, our method allows for this possibility in the form of estimating to be zero and non-zero. As demonstrated in our simulations, the loss of power in adding a single unnecessary parameter to describe many SNPs is small. Third, as discussed in the Introduction and Text S1, direct selection against disease is not a necessary condition for correlation between fitness and phenotype; as long as the disease related gene is under selective pressure in any of its functions, we expect a correlation. We have planned several extensions to this method. In addition to improved techniques of estimating fitness effects, we need to incorporate evidence for adaptive selection. Signatures of positive selection [47]–[49] can be used to prioritize genes for study which may have been more important in differentiating humans from our ancestors and hence contribute to modern phenotypes. We expect positively selected variants to have very different phenotype effects from neutral alleles, but it is not clear a-priori what that relationship should be or if it will be possible to reliably identify positively selected SNPs [50], [51]. Second, for mathematical and numerical convenience we have developed this method in the context of a prospective probability sample measuring a quantitative trait. Both these assumptions need to be relaxed for the setting of most resequencing projects. Disease phenotypes are frequently non-normal, binary, or censored such as time-to-event from clinical trials, requiring a generalized linear mixed model. The prospective sampling assumption will also require work to relax. Retrospective sampling such as in case-control designs and extreme-phenotype-based sampling [13], [52] is well known to distort random effect distributions [53]. Third, in our example and simulations, we assume that are independent of one another, but one need not do this. One could add spatial covariance structures between to relax the independence assumption, which would correspond to allowing that variants nearby each other in the genome or folded protein tend to have similar effects. Especially in exome-only resequencing studies, consideration of unobserved linked markers with techniques similar to TreeLD [35] will be important. Our model has not included dominance or epistasis between SNPs or genes, the structure of which is probably not simple, although progress has been made on determining the impact of these features to quantitative traits [54], [55]. Finally, because our example dataset comes from high-quality Sanger sequencing, we have ignored nonrandom missing data issues. Future work involving second generation sequencing or beyond must address the complex nature of library coverage, alignment error, and genotyping error inherent in those technologies.
Studies correlating genetic variation to disease and other human traits have examined mostly common mutations, partly because of technological restrictions. However, recent advances have resulted in dramatically declining costs of obtaining genomic sequence data, which provides the opportunity to detect rare genetic variation. Existing methods of analysis designed for an earlier era of technology are not optimal for discovering links to rare mutations. We take advantage of 1) the advanced theoretical understanding of evolutionary mechanics and 2) genome-wide evidence about evolutionary forces on the human genome to suggest a framework for understanding observed correlations between rare genetic variation and modern traits. The model leads to a powerful test for genetic association and to an improved interpretation of results. We demonstrate the new method on previously confirmed results in a gene related to high blood cholesterol levels.
computational biology/population genetics
genetics and genomics/complex traits
computational biology/evolutionary modeling
genetics and genomics/genetics of disease
public health and epidemiology/epidemiology
computational biology/genomics
mathematics/statistics
genetics and genomics
genetics and genomics/population genetics
An Evolutionary Framework for Association Testing in Resequencing Studies