Complex diseases and traits are likely influenced by many common and rare genetic variants and environmental factors. Detecting disease susceptibility variants is a challenging task, especially when their frequencies are low and/or their effects are small or moderate. We propose here a comprehensive hierarchical generalized linear model framework for simultaneously analyzing multiple groups of rare and common variants and relevant covariates. The proposed hierarchical generalized linear models introduce a group effect and a genetic score (i. e. , a linear combination of main-effect predictors for genetic variants) for each group of variants, and jointly they estimate the group effects and the weights of the genetic scores. This framework includes various previous methods as special cases, and it can effectively deal with both risk and protective variants in a group and can simultaneously estimate the cumulative contribution of multiple variants and their relative importance. Our computational strategy is based on extending the standard procedure for fitting generalized linear models in the statistical software R to the proposed hierarchical models, leading to the development of stable and flexible tools. The methods are illustrated with sequence data in gene ANGPTL4 from the Dallas Heart Study. The performance of the proposed procedures is further assessed via simulation studies. The methods are implemented in a freely available R package BhGLM (http: //www. ssg. uab. edu/bhglm/).
Many common human diseases and complex traits are highly heritable and are believed to be influenced by multiple genetic and environmental factors. Genome-wide association studies (GWAS) represent a powerful way for discovering disease-associated factors and investigating the genetic architecture of complex diseases [1]. In the past few years, these studies have identified hundreds of common variants (i. e. , genetic variants with minor allele frequency (MAF) >∼5%) associated with complex diseases [2]. However, the estimated effect sizes for the identified variants are small (most odds ratios are below 1. 5) and explain only a small proportion of the heritability of complex diseases [2], [3], motivating research interest in finding ‘missing’ genetic factors that contribute to the remaining heritability [4], [5]. Many explanations for the missing heritability have been suggested [4], [5]; one is that many common variants with much smaller effects are yet to be detected, and another is the possible contribution of rare variants (MAF <0. 5% or 1%) that are poorly captured by previous GWA genotyping arrays. Empirical studies and population genetics theory support the potentially important role of both rare variants and common variants of very small effects [6]–[10]. Several current studies have implicated association of rare variants with complex diseases and traits [11]–[18]. Next-generation sequencing technologies have provided unparalleled tools to sequence a large number of individuals in candidate genes, exomes, or even the entire genome, allowing for comprehensive studies of both common and rare variants. In addition to the common problems of handling large numbers of variants, however, detecting disease-associated rare variants and common variants of small effects poses unique statistical challenges [19], [20]. As such variants individually contain little variation, statistical methods that detect association between a single variant and disease phenotype provide low power with realistic sample sizes. Therefore, it is necessary to develop sophisticated methods that can effectively combine information across variants and assess the collective effect of multiple variants [4]. Several approaches along this line have been proposed [19], [20]. The basic procedure of these methods is to construct a linear combination of multiple variants with fixed weights to summarize the information across the variants and then estimate its association with the phenotype [21]–[25]. Different weights yield different summaries of the variants and implicate different assumptions about the relative importance of individual variants [24], [26]. Further, they implicitly assume that all variants affect phenotype in the same direction. However, there are many examples in which numerous rare variants detected in a gene or region may have disparate or even opposite effects on phenotype [4], [11]. Thus, these methods can be suboptimal if the data do not follow the underlying assumptions. Recently, several methods have been proposed to deal with variants with opposite effects [26]–[32], and to summarize the information across variants using non-linear functions [33], [34]. All the existing methods have been developed to assess only one group of variants at a time. Since common diseases are likely caused by a complex interplay among many genes and environmental factors, however, it is more appropriate to simultaneously model multiple groups of variants and covariates [19]. The joint analyses would improve the power of detecting causal effects and hence lead to increased understanding about the genetic architecture of diseases. Such methods are also advantageous for studies involving only one candidate gene, because numerous variants detected within a gene can be divided into multiple groups based on their allelic frequencies (common or rare) and functional annotations of the genomic regions they reside in (for example, non-synonymous or synonymous). It has been found in GWAS that the vast majority (80%) of associated variants fall outside coding regions, emphasizing the importance of including both coding and non-coding regions in the search for disease-associated variants [2]. We propose here a comprehensive hierarchical generalized linear model (GLM) framework for simultaneously analyzing multiple groups of rare and common variants and relevant covariates. The proposed hierarchical GLMs introduce a group effect and a genetic score (i. e. , a linear combination of main-effect predictors for genetic variants) for each group of variants, and jointly estimate the group effects and the weights of the genetic scores. This framework includes various previous methods as special cases, and can effectively deal with both risk and protective variants in a group and can simultaneously estimate the cumulative contribution of multiple variants and their relative importance. The methods are illustrated with sequence data in gene ANGPTL4 from the Dallas Heart Study, and are further assessed via simulation studies. Finally, we conclude this article by highlighting some areas of future research.
Suppose that a population-based association study consists of n unrelated individuals, phenotyped for a continuous or discrete disease trait and genotyped for a number of rare and/or common genetic variants in one or multiple candidate genes or genomic regions. The observed values of the response variable are denoted by y = (y1, ···, yn). We assume that the genetic variants can be divided into K groups, Gk, k = 1, ···, K, and the k-th group Gk contains Jk variants, where K≥1 and Jk>1. The groups can be constructed based on candidate genes in which the variants are located and the types of the variants (e. g. , common variants, rare non-synonymous or synonymous coding variants). We assume that some non-genetic variables (e. g. , gender indicator, age, etc.) are also measured for each individual and will be included as covariates in the model to control for possible confounding effects. We extend the hierarchical generalized linear model (GLM) of Yi and Zhi [26] to simultaneously fit covariates and multiple groups of rare and common variants. A generalized linear model consists of three components: the linear predictor η, the link function h, and the data distribution p [35], [36]. The linear predictor of individual i is expressed as the multiplicative form: (1) where β0 is the intercept, xij and βj represent covariate j and its coefficient, respectively, zij is the main-effect predictor for individual i at genetic variant j in group Gk, equaling to the number of minor alleles for an additive coding (for a rare variant, the additive coding is approximately equivalent to a dominant coding because the frequency of the minor allele is very low), the common coefficient represents the group effect for Jk variants in the k-th group, and the individual coefficients can be interpreted as the weights or relative effects of individual variants. The common coefficient represents the association between the phenotype and the linear combination of Jk individual main-effect predictors for variants in group Gk. The linear combination provides a way to combine the genetic variation across the Jk individual variants, referred to as genetic score. Therefore, the common coefficient represents the cumulative importance of the Jk individual variants in the k-th group, hence referred to as the group effect, and the weights, , give the relative importance of the individual variants in group Gk. The mean of the response variable is related to the linear predictor via a link function h: (2) The data distribution (likelihood) is expressed as (3) where is a dispersion (or variance) parameter, and the distribution can take various forms, including Normal, Gamma, Binomial, and Poisson distributions. Our main goal is to estimate the group effects and to test the hypotheses gk = 0, k = 1, ···, K. We treat the weights' s as unknown parameters and estimate them along with the group effects and other parameters from the data. But we cannot simply use classical framework (equivalent to setting uniform distributions on the' s from a Bayesian perspective), since this would result in a nonidentifiable model [37], [38]. An approach to overcoming the problem is to use an informative prior for. We use the following hierarchical prior distribution: (4) where the prior means are prefixed and will be discussed in detail later, and the subscript k[j] indexes the group k that variant j belongs to. The above hierarchical prior assumes that follows a scale mixture of normals with unknown variable-specific variance. The prior distribution for is a hierarchical formulation of the half-Cauchy distribution, which has desirable properties, such as an infinite spike at the prior mean and very heavy tails, and also facilitates efficient computation [39], [40]. An attractive feature of our hierarchical prior is that it is free of user-chosen tuning parameters and introduces group-specific parameters and variable-specific parameters and. The group-specific parameters provide a way to pool the information among variables within a group and also to induce different shrinkage for different groups, while the variable-specific parameters allow different shrinkage for different variables. Yi and Zhi [26] set the scale parameters to a known value for all the weight parameters and recommended = 0. 5 as default. However, it would be more reasonable to estimate the scale parameters from the data. If the number of groups is not large, the group effects usually can be estimated classically. However, low allelic frequencies can yield very small variances for the predictors of, i. e. , , and as a result the classical procedure can result in numerically instable estimates for the group effects. To overcome this problem, we can place a weakly informative prior on that constrains to a reasonable range [41]. We use the following hierarchical prior distribution: (5) This hierarchical prior distribution includes group-specific parameters, which can induce different shrinkage for different group effects. The group-specific parameters are assumed to follow a weakly informative prior Gamma (0. 5,0. 5). This weakly informative prior does not strongly shrink towards zero, but can constrain to lie in a reasonable range [41]. For the covariate effects, we also use the above weakly informative prior (5), i. e. , . For the intercept and the dispersion parameter, we can use any reasonable non-informative prior distributions; for example, with set to a large value, and. Our hierarchical GLMs include multiplicative parameters, a common coefficient for a group of variants and a weight parameter for each variant. As explained earlier, the common coefficient represents the overall association of the Jk individual variants in group k with the disease. In our hierarchical model, the multiplicative term can be expressed as, and thus the predictor zij ultimately gets the coefficient, which represents the main effect of that variant. The coefficient is affected by the prior mean of. Therefore, we define the adjusted main effects as, which represent the effects of individual variants. For the multiplicative model to be useful, we need informative prior distributions on the multiplicative parameters that allow us to distinguish between the group effects and the individual weights. The prior means and the variances in the normal prior distributions of the weights (i. e. ,) are the key components to interpret our hierarchical model. The variances directly control the amount of shrinkage for. If = 0, the coefficient equals the prior mean. If = ∞, is actually estimated using least squares and the parameters and cannot be distinguished. If is finite, the coefficient is shrunk towards but not identical to the prior mean. Therefore, the prior distributions bridge the gap between the two extremes of simply using the fixed weighted sum of the Jk variants as a predictor (= 0), and including them as Jk independent predictors (= ∞) [37], [38]. This interpretation can be more explicitly understood by the identity (6) where. The second term in the right side is controlled by the variances, and represents the deviation from the fixed weighted sum. Most of existing methods for analyzing rare variants proceed to construct a linear combination (genetic score) of rare variants with fixed weights [21]–[25], and thus can be viewed as special cases of our model. The prior means represent the prior relative importance of the individual variants and can be specified in various ways. The weights proposed by previous methods [21]–[25] can be used as the prior means in our hierarchical model. The simplest way is to set all = 1, resulting in the simple sum, and incorporating no prior information about the relative importance of rare variants into the model. But our method can estimate the weights from data and produce different weights to different variants based on their contributions to the phenotype. Therefore, our model uses a previous score (i. e. ,) as the baseline, and improves the fit by accounting for the variation among individual variants. An alternative choice of the prior means is to set all = 0. With this choice, the weights are shrunk towards zero, and variants with small effects can be essentially removed from the model. This seems to be reasonable, especially for the situations with non-functional variants. However, we don' t recommend this approach for rare variants for several reasons. First, most of rare and common variants have small effects, but they can be cumulatively important. In order to detect the cumulative effect, therefore, it would be better to include all the small effects in the model. Second, the estimated group effect can be less interpretable and accurate, if only one or a few variants are included in the model. Third, our hierarchical model can estimate the weights of individual variants from the data, and thus can deal with non-functional variants and disparate effects. Our Bayesian hierarchical GLMs can be fitted using Markov chain Monte Carlo (MCMC) algorithms that fully explore the joint posterior distribution by alternately sampling each parameter from its conditional posterior distribution [36]. However, it is desirable to have a faster computation that provides a point estimate (i. e. , the posterior mode) of the coefficients and their standard errors (and thus the p-values). Such an approximate calculation has been routinely applied in statistical practice [41]. We develop our mode-finding algorithm by modifying the standard iterative weighted least squares (IWLS) for fitting classical generalized linear models [42], [43]. Our algorithm updates the coefficients and using an iterative procedure. Conditional on the current estimates, we update by running the generalized linear model with the proposed prior distributions for and other corresponding parameters: (7) where, and then conditional on the current estimates, we update by running the generalized linear model with the proposed prior distributions for and other corresponding parameters: (8) where. We fit these two hierarchical generalized linear models by incorporating a flexible expectation-maximization (EM) algorithm into the iteratively weighted least squares (IWLS) for fitting classical generalized linear models. We describe our EM-IWLS algorithm in detail in Text S1. We initialize our iterative algorithm by setting the parameters () with some plausible values. For example, we can start with = 0, =, = 1, = 1, = 0. 5, = = 0. 5, and = = = = 1. We then update the parameters by iteratively running the hierarchical generalized linear models (7) and (8) until convergence. Instead of doing a nested converged EM-IWLS for each of the two models, we can run one step of the EM-IWLS at each iteration, thus taking less computing time to ultimately achieve convergence by not wasting time running many steps of the EM-IWLS within each iteration. To assess convergence, we use the standard criterion for analysis of classical generalized linear models (as implemented in the R function glm), i. e. , , where is the estimate of deviance (i. e. ,) at the tth iteration, and is a small value (say 10−5). At convergence of the algorithm, we summarize the inferences using the latest estimates of the coefficients and their standard errors. Based on these outputs, we can calculate approximate p-values as in the classical framework for testing whether a coefficient is significantly different from zero, for example, the hypothesis = 0. The adjusted main effects are then estimated as, and the approximate standard error for can be obtained by using the delta technique: (9) Therefore, we can calculate the approximate p-value for testing the hypothesis = 0. Our model fitting strategy is based on extending the well-developed IWLS algorithm for fitting classical GLMs to our Bayesian hierarchical GLMs. The IWLS algorithm is executed in the glm function in R (http: //www. r-project. org/). We have implemented the EM-IWLS algorithm by inserting the E-step for updating the missing values (i. e. , the variances and the hyperparameters and) and the steps for calculating the augmented data and the dispersion parameter into the IWLS procedure (see Text S1). We have created a new R function bglm by modifying the glm function to implement our EM-IWLS algorithm that estimates posterior modes and standard deviations for hierarchical GLMs with the prior distributions proposed here (see Text S1) and some other hierarchical priors [44], [45]. We have also developed an R function bglm. ex that implements the iterative algorithm described above for fitting our hierarchical multiplicative GLMs. Although described in the context of genetic variants in this paper, the functions bglm and bglm. ex can be used as general tools for routine data analysis using hierarchical GLMs. We have incorporated the functions bglm and bglm. ex into the freely available R package BhGLM (http: //www. ssg. uab. edu/bhglm/) that is an extensible and interactive environment for genetic association analysis of common and rare variants and gene-gene and gene-environment interactions. Our hierarchical multiplicative GLMs include various models as special cases. Although less comprehensive, these reduced models can be useful in some situations, and thus can be used as alternative approaches to analysis of multiple groups of rare and common variants. We here consider two types of reduced models. The first ignores the group effects and directly models the main effects of individual variants. Thus, the linear predictor (1) is reduced to (10) and the mean and the distribution of the response variable take the same form of the expressions (2) and (3). In this model, the coefficient represents the main effect of genetic variant j, and follows the hierarchical prior distribution (4) with the prior mean = 0. This approach can only detect individual variants with strong effects, and is less powerful in situations where the effects of all individual variants are small but they are cumulatively significant. The second alternative approach is to preset the weights of individual variants using the previous methods [21]–[25]. Thus, the linear predictor (1) becomes (11) where with fixed weights. This model is equivalent to setting the priors as (i. e. ,) and thus is a special case of our hierarchical model. The performance of this method heavily depends on the quality of the fixed weights.
We have proposed here a Bayesian hierarchical generalized linear model framework for simultaneously analyzing multiple groups of rare and common variants and relevant covariates. Since complex diseases and traits are likely influenced by multiple genetic variants and environmental factors, the joint analyses of multiple groups of genetic variants can improve the power of detecting causal effects and lead to increased understanding about the genetic architecture of diseases. The proposed hierarchical generalized linear models introduce a group effect and a genetic score for each group of variants, and jointly estimate the group effects and the weights of the genetic scores. This can produce ‘optimal’ weights to different variants based on their contributions to the phenotype, yielding an effective summary of the information across variants. The simulation studies show that the proposed method can consistently provide reasonable power even in the presence of both risk and protective variants in a group, and has better power than existing approaches even when all variants act in the same direction. Application of the method to a large published dataset on resequencing of the gene ANGPTL4 and triglycerides not only confirmed the original findings but also detected new associations. In addition to the properties described above, our method has several remarkable features. First, the proposed method can simultaneously estimate the group effects of multiple groups of variants and the individual effects of the variants, allowing us to not only identify significant genes (or groups of variants) but also assess the relative importance of single variants. Second, our hierarchical model includes various existing methods for rare variants as special cases. This shows that the proposed method is theoretically more advantageous than the existing methods, and allows us to conveniently analyze data using different ways. Third, any external information about variants, for example, the functional prediction, can be easily incorporated into our hierarchical model by specifying the prior means of the weights for variants. By doing so, our approach has the additional advantage of accounting for uncertainties about the prior assumptions. Fourth, our approach is based on the generalized linear model framework and thus can deal with various types of continuous and discrete phenotypes and covariates, and can fit any generalized linear models. Finally, the proposed algorithm extends the standard procedure for fitting classical generalized linear models in the general statistical package R to our Bayesian model, leading to the development of stable and flexible software. Our approach is highly extensible; we have planned several extensions to the proposed method, some of which have been initially implemented in our software BhGLM. The key to our approach is the use of hierarchical prior distributions for the weights and the group effects, so that these multiplicative parameters are identifiable and can be simultaneously estimated from the data. We have proposed to use the hierarchical expression of the half-Cauchy distribution with the innovation of introducing both group- and variable-specific parameters. The half-Cauchy prior is an excellent default choice for many problems [39], [40], and has been shown to perform well for our purposes. However, other hierarchical priors or penalized likelihood methods have been developed for high-dimensional data analysis, including lasso [48], [49], adaptive Lasso [50], and the elastic net [51]. These methods can be expressed as hierarchical models by assigning certain priors on the variances and other hyperparameters [45], [48], [52], and can be incorporated into our framework. Although demonstrated with only several groups of variants, our method can be adapted to deal with large-scale sequencing data involving thousands of exomes or candidate genes. For these high-dimensional settings, we need to modify the prior distributions of the group effects and the computational algorithm. We can place a shrinkage prior on the inverse scale in the gamma prior of and estimate the inverse scale from the data. We can further group the group effects based on pathways that candidate genes belong to, and specify the shrinkage priors by incorporating the second-level hierarchical structure, similar to the hierarchical priors of the weights. We describe our algorithm by simultaneously estimating all weights. This method can be very fast when the number of variables is not very large (say <2000) and has the advantage of accommodating the correlations among all the variables. However, it can be slow or even cannot be implemented when the number of variables is large due to memory storage and convergence problems. We can extend the algorithm to update coefficients group by group; at each of the iteration, the group-at-time algorithm proceeds by cycling through all the groups of parameters and treats the linear predictor of all other groups as an offset in the model. This method updates coefficients in a conditional manner, significantly reducing the number of parameters in each M-step of the EM-IWLS algorithm, and thus can deal with large number of variables. Our third extension could incorporate external gene or pathway level information into the hierarchical model. Candidate genes or pathways studies usually consist of data at different levels, i. e. , genetic variants within multiple candidate genes or pathways which may be functionally related [53]. Most of statistical methods for association studies consider only individual-level predictors (i. e. , SNPs and covariates) and ignore the hierarchical structure of the data and gene or pathway-level information. Often, rich gene or pathway-level information is available [54], including genomic annotation or pathway ontologies, functional assays, in silico predictions of function or evolutionary conservation [55]. Therefore, there is a growing need to develop sophisticated approaches that model the multilevel variation simultaneously and incorporate gene or pathway-level data into the model [56], [57]. Our hierarchical models provide a natural and efficient way to incorporate the external information about candidate genes into the analysis. One way to include the gene-level information in the hierarchical models is to model the prior means of weights and group effects using gene or pathway level predictors [38]. This would allow us to pool the information in the same genes or pathways and thus would provide more effective inference about the genetic effects. Our fourth extension could incorporate genetic interactions (gene-gene and gene-environment interactions) into the model. Just as interactions must be considered in standard GWA studies [57]–[59], they are also likely to be important in association studies involving rare variants [19]. In principle, we can extend the proposed model to include additional groups for interactions for each pair of groups of main effects and to define an overall effect and a genetic score for each interaction group. However, it would be necessary to investigate statistical power for detecting interactions for rare variants. Finally, we have planned to extend our method to family-based matched case-control association studies. So far the existing methods for rare variants have focused on population-based studies. However, for rare variants, family-based designs may prove very useful [60]. Not only are they robust against population stratification, but they may also offer increased power due to the increased likelihood of affected relatives to share the same rare disease variants. As the conditional logistic regression commonly used for matched case-control studies can be formulated as a Poisson regression [36], our hierarchical generalized linear models can be applied. The proposed hierarchical generalized linear models may provide efficient tools for disease risk prediction and personalized medicine. GWA studies have raised expectations for predicting individual susceptibility to common diseases using genetic variants [61], [62]. Previous methods using only a limited number of significant variants have typically failed to achieve satisfactory prediction performance [63], [64]. Recent studies show that joint analysis of a large number of genetic variants can improve the prediction of complex traits [65]–[67]. It is understood that a model including as many predictors as possible and fitted appropriately could provide better prediction. Although the previous studies have included many genetic variants in a predictive model, they treat these variants individually and hence could be suboptimal to efficiently use information of genetic variants with small effects and low frequencies. The proposed hierarchical models can better deal with such variants and can integrate external biological knowledge, and therefore may be able to improve the accuracy of prediction.
Complex diseases and traits are likely influenced by many common and rare genetic variants and environmental factors. Next-generation sequencing technologies have provided unparalleled tools to sequence a large number of individuals, allowing for comprehensive studies of both common and rare variants. However, detecting disease-associated rare variants and common variants of small effects poses unique statistical challenges. We propose here a comprehensive hierarchical generalized linear model framework for simultaneously analyzing multiple groups of rare and common variants and relevant covariates. The proposed hierarchical generalized linear models introduce a group effect and a genetic score for each group of variants, and jointly they estimate the group effects and the weights of the genetic scores. This framework includes various previous methods as special cases, and it can effectively deal with both risk and protective variants in a group and can simultaneously estimate the cumulative contribution of multiple variants and their relative importance. The methods are illustrated with sequence data in gene ANGPTL4 from the Dallas Heart Study and are further assessed via simulation studies. The methods have been implemented in a freely available R package BhGLM (http: //www. ssg. uab. edu/bhglm/).
medicine
computer science
mathematics
genetics
biology
genetics and genomics
Hierarchical Generalized Linear Models for Multiple Groups of Rare and Common Variants: Jointly Estimating Group and Individual-Variant Effects