PDF

J Genomics 2014; 2:1-19. doi:10.7150/jgen.5054 This volume Cite

Research Paper

Meta-Analysis of Candidate Gene Effects Using Bayesian Parametric and Non-Parametric Approaches

Xiao-Lin Wu^1,2 , Daniel Gianola^1,2,3, Guilherme J. M. Rosa^2,3, Kent A. Weigel¹

1. Department of Dairy Science, University of Wisconsin, Madison, WI 53706, USA;
2. Department of Animal Sciences, University of Wisconsin, Madison, WI 53706, USA;
3. Department of Biostatistics and Medical Informatics, University of Wisconsin, Madison, WI 53706, USA.

✉ Corresponding author: Xiao-Lin Wu, 1675 Observatory Dr., Department of Dairy Science, University of Wisconsin, Madison, WI 53706, USA. Tel.: (608) 263 7824; Fax: (608) 263-9412; E-mail: nick.wuwisc.edu.

Abstract

Candidate gene (CG) approaches provide a strategy for identification and characterization of major genes underlying complex phenotypes such as production traits and susceptibility to diseases, but the conclusions tend to be inconsistent across individual studies. Meta-analysis approaches can deal with these situations, e.g., by pooling effect-size estimates or combining P values from multiple studies. In this paper, we evaluated the performance of two types of statistical models, parametric and non-parametric, for meta-analysis of CG effects using simulated data. Both models estimated a “central” effect size while taking into account heterogeneity over individual studies. The empirical distribution of study-specific CG effects was multi-modal. The parametric model assumed a normal distribution for the study-specific CG effects whereas the non-parametric model relaxed this assumption by posing a more general distribution with a Dirichlet process prior (DPP). Results indicated that the meta-analysis approaches could reduce false positive or false negative rates by pooling strengths from multiple studies, as compared to individual studies. In addition, the non-parametric, DPP model captured the variation of the “data” better than its parametric counterpart.

Keywords: Bayesian models, candidate genes, Dirichlet process prior, Markov chain Monte Carlo, meta-analysis.

Introduction

Candidate gene (CG) approaches constitute a complementary strategy to map-based cloning and insertional mutagenesis [1-2]. A candidate gene is either a cloned gene presumed to be involved in the regulation and expression of a given trait (“functional candidate”) or a gene closely linked to loci controlling the trait (“positional candidate”). An appealing feature of the CG approach is that it does not require the development of mapping populations, and the analysis relies simply on association tests between polymorphisms of a molecular marker (or markers) and variation for the trait of interest.

Candidate gene approaches have been widely used for identification and characterization of Mendelian and quantitative trait loci (QTLs) in the past decade [3], because they are simple and convenient. The number of published candidate gene studies has been increasing [4]; In pigs, for example, a search in PubMed (http://www.ncbi.nlm.nih.gov/sites/entrez) shows over 40 published papers about candidate genes on litter size, since the work of Rothschild et al. (1996) [2]. This number represents those deposited in the PubMed databases, which makes up only a portion of the studies that have been conducted so far in search for genes affecting pig litter size.

The large number of candidate gene studies, however, does not provide an indication of the reliability of these results, and their conclusions tend to be inconsistent among individual studies. Frequently, when a candidate gene is studied by several groups or even by the same group using different experimental populations, each association found previously can be paired with an equally convincing rebuttal, and both positive and negative findings are likely to be observed in reality. Even when a genuine genetic association exists, the effect size of the gene varies among different studies, hence leading to various conclusions. Such inconsistencies may result from differences in many issues, which include and not limited to varied populations, methodologies, and experimental designs used in CG studies.

This aforementioned situation can be handled by meta-analyses, which typically involves pooling effect-size estimates or combining P values from multiple studies to arrive at an overall conclusion [5,6]. Both fixed-effects and random-effects models have been proposed for meta-analyses [3]. Often, the model poses a center-specific “true effect” while taking into account the heterogeneity over individual studies [7,8]. Hence, meta-analyses provides a powerful tool for assessing population-wide effects of candidate genes on quantitative traits and it can reveal diversity or heterogeneity that are previously undiscovered [9,10].

Current meta-analytic random-effects models assume that the CG effect size among studies is normally distributed, which may not be true in reality. For instance, in a meta-analysis evaluating the association of the platelet PIA polymorphism of glycoprotein IIIa and risk of coronary heart disease, Burr et al. (2003) found that kernel density estimates of study-specific CG effects were multi-modal [7], which obviously violated the normality assumption. In reality, population stratification and admixture is common, leading to multi-modal distributions of the data [11].

In the present paper, we first show how a mixed-effect model can be used to describe the effects of a candidate gene obtained from multiple studies, and then we extend this meta-analytic model to non-parametric settings. The parametric model assumes a normal distribution of study-specific gene effects (referred to as the meta-NP model) whereas the non-parametric model relaxes the normality assumption by replacing it with a more general, discrete, “nonparametric prior” - the Dirichlet process prior (referred to as the DPP model). Both models are implemented within the Bayesian framework. Features of these methods are illustrated using simulated data.

Statistical Methods

Mixed models for study-specific CG effects

Consider N studies that investigate the effect of a candidate gene on a quantitative trait (the methodology can be extended to a set of candidate genes evaluated jointly). Let y_i be a point estimate of the effect size u_i of the candidate gene (e.g., additive or dominance effects, measured as a linear function of genotypic values) in the i-th study, each representing a single population. Then, y_i's can be modeled as:

y_i | u_i, σ_i² ~ N (u_i, σ_i²), for i=1,…,N.  (1) 

Where σ_i² is assumed to be known but empirically equal to the sample variance, s_i² in the i-th population. The assumption σ_i² = s_i² is typically made in meta-analysis, which ensures the parameter identification in the model. This assumption is valid as long as sample sizes in individual studies are not small [7,12].

Specifically, if we let y_i be the estimated additive effect of the candidate gene in the i-th population, computed as one half of the difference between mean values of the two homozygotes, then, assuming no interactions between this gene and other genes also affecting the trait of interest, we have

 (2) 

where, for examples, is an estimate of the genetic effect of genotype AA, pertaining to the i-th study, a_AA is its true genetic value, is the average residual genetic values of individuals having genotype AA, and is the average of the errors for all individuals having genotype AA at the CG locus. The residual genetic values are attributable to loci other than the candidate gene under investigation. The assumption made in (2) is that the candidate gene is the functional gene affecting the trait. If the candidate gene is only a marker linked to a functional gene, its association (additive) effect, denoted as a*, can be expressed as a function of the gene effect a and of the recombination rate r between the marker and the causal gene, e.g., a* = a(1-2r).

Evidently, if and , then provides an unbiased estimator of (a_AA - a_aa). This situation, however, is unlikely to be true because it rests on over-simplified assumptions. For instance, the genetic determination of a quantitative trait may involve some major genes as well as polygenic effects, and occasional associations between these genes and the quantitative trait can lead to misinterpretation of estimated candidate gene effect, unless they are accounted for properly. Assuming that the quantitative trait is additively affected by J major genes, as well as infinitesimal gene effects, then, without considering interactions among them, can be decomposed as:

 (3) 

where, for example, b_AA,j is the effect of genotype AA for the j-th major gene, f_AA,j|AA(i) is the frequency that the genotype of the j-th major gene is AA for all individuals whose genotype at the CG locus is also AA, and is the average value of polygenes of individuals having genotype AA in the i-th population (candidate gene study). Then, if we denote and , the first term b_i^* reflects genetic “noise” due to some major genes affecting the trait, whereas ω_i^* is attributable to polygenic effects. Thus, it becomes apparent that assuming may not be as reasonable as one can expect. In reality, due to effects of major genes and polygenes that are not randomized among the three CG genotype groups.

Now, let's revisit (2). Denote a_AA-a_aa=2μ, which is invariant across individuals populations. Thus, we have:

y_i = μ + g_i + e_i  (4) 

where . Here, we assume e_i ~ N(0, s_i²), where s_i² is approximated by the variance of the estimate in study i. This decomposition of candidate gene effects suggests that a mixed-effect model can be used to model candidate gene effects estimated from multiple studies. In (4), μ corresponds to the “central” putative effect of the candidate gene, and g_i represents random study-specific deviations from the “central” effect size in the i-th study with some variance. Therefore, as long as the genetic background is not randomized properly in a CG study, we have E(g_i) ≠ 0 and estimated CG effect in an individual CG study is biased even if this candidate gene represents the true function mutation of the QTL.

In a meta-analysis, however, the estimation of the “central” effects involves weighting estimates from single studies across a variety of genetic backgrounds. This can be viewed as kind of randomization of the genetic background effects across multiple experiments in some sense, which in turn can minimize the bias in estimated CG effects. Re-write (4) in matrix notation as:

 (5) 

Where y = (y₁ y₂ … y_N)', g = (g₁ g₂ … g_N)', and is an N×N known (diagonal) covariance matrix with being the i-th diagonal element with all off-diagonal elements being equal to zero. Further, we assume that the g_i's follow the same normal distribution with mean zero and variance σ_g². That is,

 (6) 

Here, σ_g² is the variance of study-specific deviations of the CG effect from the overall mean μ due to background major genes and polygenes.

To complete the Bayesian settings, a normal and a scaled inverse chi-square prior distribution are assumed for μ and σ², respectively. That is, p(μ) = N (μ₀, σ_μ²) and , where μ₀, σ_μ², υ_g and Ѕ_g² are known hyper-parameters.

Denote θ = (μ, g, σ_g²), and let H represents all known hyperparamters. Then, the joint distribution of θ given the data y and hyperparameters H is

 (7) 

The joint posterior distribution can be evaluated by using a Markov chain Monte Carlo approach such as Gibbs sampling where posterior samples are drawn from the fully conditional distributions of the model parameters, either element-by-element or block-wise [13]. Briefly, the conditional posterior distribution of μ is the normal distribution

 (8) 

where:

and

In (8), “else” is used to represent all parameters other than μ. The conditional posterior distribution of g is multivariate normal:

 (9) 

where and . Finally, the fully conditional distribution of the genetic variance σ_g² is the scaled inverse Chi-squared distribution:

 (10) 

Bayesian non-parametric hierarchical model

There are essentially two normality assumptions in the parametric model. The first one, as seen in (1), relies on the central limit theory. In other words, the sampling distribution of the estimate of a candidate gene effect tends to be a normal distribution as the sample size increases. However, the second normality assumption about the distribution of study-specific CG effects, as in (4), does not have any theoretical justification. For example, Burr et al. (2003) gave evidences of non-normality of CG effects among a number of independent studies [7].

Considering that the candidate gene effect sizes over multiple studies can departure from a homogeneous normal distribution, we replace (6) with , , where G is some general distribution, G ~ π, and π is a “nonparametric prior”. A choice for π is the Dirichlet process [14-16]. That is,

G ~ DP(αG₀)  (11) 

Above, G₀ is a “baseline” distribution function, which is also referred to as the “center” of the DP prior, in the sense that for any given θ we have E(G(θ)) = G₀ (θ), and α is interpretable as a precision parameter indicating the degree of concentration of the prior on G around some parametric family { G₀(θ); θ∈Θ}. We assume G₀ (σ_g²) ≡ N (g_i |0, σ_g²), where σ_g² is the variance of the baseline distribution G₀.

Denote θ ={ μ, g, σ_g²}, and let H contain all known hyper-parameters. In the first stage, it is assumed that the value of α is known. So, the joint conditional distribution of θ, given y, H, and is

 (12) 

Above, g ~ G denotes that the g_i's follow the unknown distribution G, drawn from a Dirichlet process with G₀ ≡ N (g_i |0, σ_g²) as base measure and known precision parameter μ. The notation [g_k ~ G][G ~ DP(αG₀(μ, σ_g²))] indicates that the distribution G is integrated out using the DP process as mixing entity. A priori, μ and σ_g² follow the same distributions as their counterparts in the parametric model discussed previously.

Bayesian implementation of this meta-DPP model via Markov chain Monte Carlo can be conducted following Gianola et al. (2010) [17] and Wu et al. (2011) [18]. Briefly, given everything else, the posterior distribution of μ in the non-parametric model is normal with the same form as in (8). Since the posterior distribution of μ is not changed when switching to a DP prior, sampling of μ is as in the parametric model. But sampling the g_i's in the non-parametric model is different from that in the parametric model. The conditional posterior distribution of g_i, given anything else, can be expressed as

 (13) 

Here, because random distribution G has been integrated out at this point, one does not need to sample G. In (13), we observe that that the fully conditional posterior distribution of g_i is a mixture of N-1 degenerate distributions, δ(g_i'), with point mass at g_i, (i'≠i), and of the parametric conditional distribution under the assumption of a baseline distribution. Specifically, if α→∞, with G₀ ≡ N(0, σ_g²), then (13) tends to (9). In other words, the posterior distribution of g in the non-parametric model is approaching a normal distribution as α→∞.

In essence, the posterior distribution of g is discrete, implying a clustering property of this meta-DPP model. This in turn groups candidate gene effects with similar effects from multiple studies. Let γ₁,…, γ_k, be K distinct values among N estimates of the gene effect size g₁,…, g_N, where 1≤K≤N, at the t-th iteration of the sampler. Then, the unique values of g₁,…, g_N induce a partitioning into K clusters. For each cluster, say k, the g_j's take on the same value γ_k. Given that g₁,…, g_N are random, this induces a random portion as well [17,18]. In the Bayesian context, the marginal posterior distribution of the number of clusters can be constructed, e.g., by directly counting the number of distinct values of g_i's at each iteration of the sampler. Following Bush and MacEachern (1996), we re-sample the γ's after the groups of effects are determined, which improves mixing [19].

Assuming a scaled inverse chi-square prior distribution, σ_g² ~ χ^-2(υ_g,S_g²), the conditional posterior distribution of σ_g² is also scaled inverse chi-square:

 (14) 

Finally, the precision parameter α can be treated as unknown and estimated in the analysis. Following West (1992) [20], given a Gamma prior distribution, p(α) = Gamma(a,b), the fully conditional distribution of α is a mixture of two Gamma distributions, with the mixing probabilities being π_s and 1 - π_s, as:

 (15) 

where a* = K + a, b* = b- log s, and . The auxiliary variable s is drawn from a Beta distribution (Gianola et al., in press):

p(s|K, α, else) = Beta (α+1, K)  (16) 

Simulation studies

QTL, candidate gene, and the quantitative trait

Five meta-analyses were simulated, each consisted of 30 individual studies. In principle, each study represented a specific population under investigation, with the number of individuals per population (sample size) ranging from 52 to 498 (Figure 1). Among them, four populations consisted of less than 100 individuals (i.e., 71, 57, 52, 53 individuals, respectively), mimicking real situations in which samples in large numbers can be difficult to obtain.

 Figure 1 

Number of individuals per population in the simulation studies.

The simulation studies were conducted under two scenarios. The first scenario, as in the first meta-analysis, assumed that the quantitative trait was affected by one bi-allelic QTL and infinitesimal gene effects as well. The additive effect for the QTL was 0.8 units, and the dominance effect of the QTL was non-existent. Essentially, this was an additive QTL. For the sake of simplicity, a bi-allelic candidate gene was simulated, which exactly coincided with the QTL. In other words, the candidate gene fully represents a functional mutation of the QTL. The effect size of the candidate gene alleles followed the order A > a. The allelic frequencies (say allele a) of the candidate gene varied in the 30 populations, ranging from 0.166 to 0.845 (Figure 2). Assuming random mating, instances of the three genotypes at the candidate gene locus were generated from the multinomial distribution

where f_i is the frequency of allele a and, for example, n_aa is the number of individuals having genotypes aa in the i-th population, with N = n_aa + n_Aa + n_AA.

The second scenario was represented by four simulation studies (denoted by S2a, S2b, S2c, and S2d), in which the quantitative trait was genetically determined by two QTLs, one additive and one dominant, and by infinitesimal genetic effects as well. The two QTLs were not linked. The additive effect of both QTLs was 0.6. The dominance effect was 0 for the first QTL and it was 0.3 for the second QTL. Again, the candidate gene fully represented a functional mutation of the first QTL affecting the quantitative trait. The allele frequencies of the first QTL were comparable among populations, which was approximately 0.4 for allele a and 0.6 for allele A, but the allele frequencies of the second QTL varied in the 30 populations, ranging between 0 and 1 (Figure 3). The genotypes of a QTL (candidate gene) were sampled from a multinomial distribution with probability, say, p(G=aa) = f_i², p(G=Aa) = 2f_i (1- f_i), and p(G=AA) = (1- f_i)², where f_i is the frequency of allele q (or a) at the i-th population. Under both scenarios, the overall mean was the only fixed effect, which was arbitrarily given. The infinitesimal genetic effects of each individual followed a normal distribution with mean zero and variance 1.0. The residuals of the phenotypic values followed a normal distribution with mean zero and variance 4.

 Figure 2 

Frequencies of allele a at the candidate gene locus in 30 populations, as generated in simulation S1.

 Figure 3 

Frequencies of allele q of the second QTL (not the candidate gene) affecting the trait in the 30 populations, as generated in simulations S2a, S2b, S2c, and S2d. The candidate gene is the first QTL itself, whose allelic frequencies are on average 0.4 for allele a and 0.6 for allele A, respectively.

The quantitative trait for each individual was simulated as the sum of the overall mean, QTL effects, infinitesimal effects, and residuals. Systematic environmental factors were not considered in this series of simulation studies, yet they can be relevant in practical situations. The data used in the subsequent meta-analysis consisted of additive and dominance effects of the candidate gene estimated from a total of 30 populations, which is described as follows. For each population, the additive effect of the candidate gene was computed as:

 (17) 

where and were means of the quantitative trait of all individuals having genotype AA and aa, respectively. Assuming homogeneous variance for the three genotypes of the candidate gene, the variance is estimated by:

 (18) 

where, for example, n_aa(i) and S²_aa(i) are the number of individuals and the sample variance of genotype aa. Then, the standard error of the additive effects is calculated as:

 (19) 

The dominance effect of the candidate gene and the standard error in the i-th study are computed, respectively, by:

 (20) 

and

 (21) 

where is the trait mean of all individuals having genotype Aa and n_Aa(i) is the corresponding sample size. In addition, additive and dominance effects of these candidate genes were also estimated using the linear regression approach. Briefly, the genotypes are coded as -0.5 (aa), 0 (Aa), and 0.5 (AA) for the additive effect, and -0.5 (aa), 1 (Aa), and -0.5 (AA) for the dominance effect. Regression coefficients were estimated using the lm function provided by the R Stats package [21].

Markov chain Monte Carlo sampling

The meta analysis was implemented by Bayesian analysis via Markov chain Monte Carlo (MCMC) sampling, based on the meta-NP and the meta-DPP models, respectively. Each analysis consisted of 200,000 iterations after a burn-in period of 2,000 iterations. The posterior samples were thinned every one-tenth and saved for posterior inference of model parameters. Convergences of the chains were diagnosed visually, and all the chains converged and mixed well.

Results and Discussion

Estimation of candidate gene effects in multiple studies

For each of the 30 CG studies, additive and dominance effects (and their standard errors) of the candidate genes were calculated using formula (17)-(21), and presented in Tables 1 and 2, respectively. These estimates corresponded well to those obtained using the linear regression approach. For example, plots of linear fitting between estimated additive and dominance effects of the candidate gene (and their standard errors) obtained using both approaches in meta-analysis S2a are illustrated in Figure 4. The R², which measures the goodness of fit between estimates obtained using both approaches, was approximately 1 for both additive and dominance effects, and for their standard errors.

Conclusions from individual CG studies varied drastically. In meta-analysis S1, for example, additive effects of the candidate gene varied dramatically among the 30 independent studies, ranging from 0.059 to 1.630 (Table 1). The differences reflected largely the consequence of genetic sampling of the QTL and infinitesimal genes. While most independent studies suggested the presence of significant (p<0.05) or very significant (p<0.001) additive effect of the candidate gene, there were several studies (i.e., from 5 to 10 studies in each of the meta-analysis) showing non-significant (p>0.05) results. In S2d, even negative estimates of additive effect for the candidate gene were obtained in three independent studies, though, by simulation, the additive effect was positive.

The simulated dominance effect of the candidate gene (corresponding to the first QTL) was zero, and estimated dominance effects (Table 2) were not significant (p>0.05) in most of the 30 independent studies. However, as shown in Table 2, there were still five studies supporting the presence of significant dominant effect, with the p values ranging from 0.010 to 0.035.

The above results clearly demonstrated that individual candidate gene studies tend to yield high rates of false positive and false negative results, leading to inconsistent conclusions among individual studies. As was mentioned previously, estimated candidate gene effects are biased if E(g_i) ≠ 0, and the bias increases with the difference of allele frequencies of other functional loci affecting the trait. In reality, the bias can result from many factors including linkage and assortative mating, also possibly from population stratification and admixture. In the present simulation studies, genetic sampling of the QTL and infinitesimal gene effect was the major cause of the bias.

 Figure 4 

Comparison of candidate gene effects and their standard errors obtained using formula (18)-(21) (represented by the y-axis) and those obtained using the linear regression approach (represented by the x-axis) in meta-analysis S2a. These graphs represent: (a) additive effects, (b) standard errors of additive effects, (c) dominance effects, and (d) standard errors of dominance effects.

Bayesian parametric meta-analysis of candidate gene effects

In the meta-NP model, in which a normal prior distribution was assumed for the study-specific CG effect, the overall mean estimated the “central” candidate gene effect and it corresponded well to the simulated value of the candidate gene. The overall mean of additive effect of the candidate gene was estimated to be 0.811 in meta-analysis S1 and it was from 0.540 to 0.631 in meta-analyses S2a, S2b, S2c and S2d (Table 3). These estimates agreed with the simulated values of the candidate gene, which were 0.8 in S1 and 0.6 in S2a-S2d, respectively. The 95% highest posterior density (HPD) interval of the additive effects obtained from the meta-NP model always contained the “true” values. Estimated overall mean of dominance effect of the candidate gene was small (for instance, -0.092 in S1 and from -0.124 to -0.021 in S2a-S2d) and not significant because the 95% HPD intervals contained zero (Table 4).

However, the normality assumption for the study-specific deviations of the candidate gene effects is questionable. As shown in Figures 5 and 6, kernel density plots of additive and dominance effects of the candidate gene (e.g., in simulations S1, S2b, and S2d) were multi-modal, and them by no means looked like “bell-shaped” normal distributions. This finding coincided with that of Bur et al (2003) [7].

Bayesian non-parametric meta-analysis of candidate gene effects

The meta-DPP model relaxed the normality assumption for study-specific CG effect by replacing it with a general prior distribution - the DP prior. Noticeably, the meta-DPP model produced estimates of the “central” candidate gene effects that were generally closer to the simulation values and with smaller standard errors (Table 3), as compared to the parametric method with a normal prior (referred to as the meta-NP method). In S2a, for example, the estimate (standard error) for the mean of the additive effect of the candidate gene was 0.569 (0.058) from the meta-DPP method and 0.540 (0.074) from the meta-NP method. Same patterns were observed in simulations S2b, S2c, and S2d (Table 3). Nevertheless, the overall means of the additive effect from both approaches were similar to each other and they both corresponded well to the “true” value. This high agreement between the two models was probably because the candidate gene was the only major gene (QTL) affecting the quantitative trait in S1, and the effect of genetic sampling was not as dramatic as it was in simulations S2a-S2d when a second QTL was also acting. Even in the latter case, the DPP estimates still had a smaller standard error than the meta-NP estimate. In all five simulations, the 95%HPD interval for the additive effect contained the true values. The posterior distributions of the “central” additive effects were approximately, normally distributed and symmetric, and examples from meta-analyses S1, S2a, and S2c are illustrated in Figure 7.

 Figure 5 

Kernel density plots of additive effects of the candidate gene obtained from 30 independent studies in meta-analyses S1, S2b, and S2d. The X-axis represents additive effect and the Y-axis represents kernel density estimates.

 Figure 6 

Kernel density plots of dominance effects of the candidate gene obtained from 30 independent studies in simulations S1, S2b, and S2d. The X-axis represents dominance effect and the Y-axis represents kernel density estimates.

 Figure 7 

Posterior distributions of “central” additive effect of the candidate gene in simulations S1, S2b, S2d, each representing a meta-analysis with 30 individual studies (i.e., 30 independent populations). The x-axis represents values of posterior draws and the y-axis represents kernel density estimates.

The meta-DPP model captured the feature of the data well because it relaxed the normal assumption about the distribution of the candidate gene effects. In Figure 8a, the histogram and kernel density plots of the estimated additive effects from 30 individual studies in meta-analysis S2c clearly shows a multi-modal distribution, rather than a normal distribution typically assumed in the meta-NP model. Figure 8b shows an instance of initial states of cluster values of study-specific effects (specifically, study-specific effects plus the “central” additive effect μ, such that they are on the same scale as the estimates from the individuals studies), and Figure 8c shows an instance of cluster values during MCMC sampling (the “central” additive effect μ was also added for the sake of comparison). Apparently, the clustering property of the meta-DPP method induced an aggregated structure of the effect sizes of the candidate gene from single studies, which could capture the variation of candidate gene effects obtained from multiple studies.

The meta-DPP approach suggested that the dominance effect of the candidate gene was not significant, because the 95% HPD interval contained zero (Figure 9). This was in agreement with the meta-NP approach. Both types of meta-analyses could reduce substantially the chance of reporting false positive or false negative candidate effects on the quantitative trait, as compared to individual candidate gene studies.

Posterior distributions of the variance of the study-specific CG effects (σ_g²), as obtained from the meta-DPP method, are illustrated in Figure 10. The estimates (standard deviations) of σ_g² obtained from the meta-DPP method were larger than those from the meta-NP method (Tables 3 and 4).

Summary posterior statistics of other model parameters (i.e., K, α, s) are listed in Table 5. On average, the meta-DPP approach postulated from 1.6 to 3.0 clusters for the additive effect, and from 1.7 to 3.5 clusters for the dominance effects in the five simulation studies. The posterior distributions of the number of clusters, e.g., for additive effects of the candidate gene in meta-analysis S1, S2b and S2d, were shown in Figure 11. The posterior means of α were small, ranging from 0.705 to 1.026 for the additive effect and from 0.739 to 1.139 for the dominance effect. This was indication that the data were diffuse, and the distribution of candidate gene effects among studies deviated from a normal distribution. The small values of α might explain the difference in the estimates of σ_g² between the two the parametric and non-parametric models. As α becomes large, the distribution of σ_g² in the meta-DPP model would approach a normal distribution, and the estimates of σ_g² would become comparable between the two approaches.

Finally, the clustering of g_i's by their inferred values is a noting feature of the non-parametric method with DP prior. In meta-analysis S2b, for example, when the additive effects were initialized as belonging to over 10 clusters, they could converge to two clusters (Figure 12a,b). Similarly, when the dominance effects (e.g., in simulation S2d) were initialized with over 10 clusters, they could converge to approximately three clusters or less (Figure 12c,d). Theoretically, the meta-DPP model can be considered as an infinite model, that is, a mixture model with a countably infinite number of clusters. In practice, however, the probability of drawing g_i is decreased quickly, and only a small number of clusters are possibly. Intuitively, imaging we generate samples from . Because G ~ DP(αG₀) is discrete with positive probability, there will be ties among the g_i's. So, the g_i's will form clumps. When α is small, the first few probabilities P_i's add up to nearly 1, resulting in higher probability ties. This leads to important consequences regarding the posterior distribution of g given the data y. Consider, for example, the distribution of g_i. Conditional on the data and g_-i (i.e., g without g_i), the estimation results in shrinkage towards and towards a grand mean. Because of the propensity of clumping, the posterior is also shrunk towards those in g_-i whose values are close to . This property results in a way of pooling information that involves weighing results of similar studies more heavily [22].

In summary, we have presented and compared two meta-analysis models that were used to evaluate candidate gene effects using simulation data. While both models applied well in the meta-analysis of candidate gene studies, the non-parametric model with DPP captured better the distribution of the candidate gene effects among studies, hence leading to more robust conclusions. The non-parametric model with DPP makes weaker assumption about the distribution of candidate genes among individual studies, and hence it model can outperform the parametric model in fitting the data with cryptic gene interaction. Finally, it should be noted that a candidate gene may not be the causative locus itself, but linked gene. To further identifying causative loci is beyond the scope of the present study. Also, with two QTLs, we did not evaluation their interactions in the present study. Instead, we retain simple QTL settings to compare the performance of the two meta-analysis models.

 Figure 8 

Histograms and kernel density plots of: (a) additive effects of the candidate gene obtained from 30 studies, (b) initial state of the cluster values, and (c) a state of the cluster values during the MCMC sampling, in simulation S2c.

 Figure 9 

Posterior distributions of “central” dominance effect of the candidate gene in meta analyses S1, S2b, and S2d, respectively, each consisting of 30 independent studies. The x-axis represents values of posterior draws and the y-axis represents kernel density estimates.

 Figure 10 

Posterior distributions of the variance of study-specific candidate gene (additive) effects in meta analyses S1, S2b, and S2d, respectively, each consisting of 30 independent studies. The x-axis represents values of posterior draws and the y-axis represents kernel density estimates.

 Figure 11 

Posterior distributions of the number of clusters of additive effects of the candidate gene in meta analyses S1, S2b, and S2d, respectively, each consisting of 30 independent studies. The x-axis represents values of posterior draws and the y-axis represents kernel density estimates.

 Figure 12 

Comparison of states of cluster values assigned to the 30 population at the initial stage (left) and a stage after convergence (right), respectively. Graphs (a)-(b) correspond to the model with additive effect in meta-analysis S2b; graphs (c)-(d) correspond to the model with dominance effects in meta-analysis S2d.

Acknowledgements

This research is supported by the Wisconsin Agriculture Experiment Station, and by grant NSF DMS-044371.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Pflieger S, Lefebvre V, Causse M. The candidate gene approach in plant genetics: A review. Mol Breed. 1996;7:275-91

2. Rothschild M, Jacobson C, Vaske D. et al. The estrogen receptor locus is associated with a major gene influencing litter size in pigs. Proc Natl Acad Sci U S A. 1996;93:201-5

3. Wu XL, Hu ZL. Meta-analysis of QTL Mapping Experiments. Methods Mol Biol. 2012;871:145-71

4. Munafò MR. Candidate gene studies in the 21st century: meta-analysis, mediation, moderation. Genes Brain Behav. 2006;5(Suppl 1):3-8

5. DerSimonian R, Laird N. Meta-analysis in clinic trials. Control Clin Trials. 1986;7:177-88

6. Normand ST. Meta-analysis: formulating, evaluating, combining, and reporting. Stat Med. 1999;18:321-59

7. Burr D, Doss H, Cooke GE. et al. A meta-analysis of studies on the association of the platelet PlA polymorphism of glycoprotein IIIa and risk of coronary heart disease. Stat Med. 2003;22:1741-60

8. Burr D. Doss H. A Bayesian semiparametric model for random-effects meta analysis. J Amer Stat Assoc. 2005;100:242-51

9. Ioannidis JP, Ntzani EE, Trikalinos TA. et al. Replication validity of genetic association studies. Nat Genet. 2001;29:306-9

10. Ioannidis JP, Trikalinos TA, Ntzani EE. et al. Genetic associations in large versus small studies: an empirical assessment. Lancet. 2003;361:567-71

11. Clayton D. Population association. In: (ed.) Balding D J, Bishop M, Cannings C. Handbook of Statistical Genetics. Chichester: John Wiley & Sons Ltd. 2001:519-40

12. DuMouchel W. Bayesian meta analysis. In: (ed.) Berry D. Statistical Methodology in the Pharmaceutical Science. New York: Marcel Dekker. 1990:509-29

13. Sorensen D, Gianola D. Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. New York: Springer-Verlag. 2002

14. Ferguson TS. A Bayesian analysis of some nonparametric problems. Ann Stat. 1973;1:209-30

15. Ferguson TS. Prior distributions on spaces of probability measures. Ann Stat. 1974;2:615-29

16. Antoniak C. Mixtures of Dirichlet process with applications to Bayesian nonparametric problems. Ann Stat. 1974;2:1152-74

17. Gianola D, Wu XL, Manfredi E. et al. A non-parametric mixture model for genome-enabled prediction of genetic value for a quantitative trait. Genetica. 2010;138:959-77

18. Wu XL, Gianola D, Hu ZL. et al. Meta-Analysis of Quantitative Trait Association and Mapping Studies using Parametric and Non-Parametric Models. J Biomet Biostat. 2011 S1

19. Bush CA, MacEachern SN. A semiparametric Bayesian model for randomized block designs. Biometrika. 1996;83:275-85

20. West M. Hyperparameter estimation in Dirichlet process mixture models. Duke University Technical Report 92-A03. 1992

21. R: A language and environment for statistical computing. R Development Core Team. http://www.R-project.org/

22. Escobar MD, West M. Bayesian density estimation and inference using mixtures. J Amer Stat Assoc. 1995;90:577-88

Author contact

Corresponding author: Xiao-Lin Wu, 1675 Observatory Dr., Department of Dairy Science, University of Wisconsin, Madison, WI 53706, USA. Tel.: (608) 263 7824; Fax: (608) 263-9412; E-mail: nick.wuwisc.edu.

Published 2014-1-2

Citation styles

©2024 Ivyspring International Publisher. Terms of use