A Hierarchical Bayesian Model for a Novel Sparse Partial ... - Genetics

Copyright Ó 2010 by the Genetics Society of America DOI: 10.1534/genetics.110.115055

A Hierarchical Bayesian Model for a Novel Sparse Partial Diallel Crossing Design Anthony J. Greenberg,*,1 Sean R. Hackett,* Lawrence G. Harshman† and Andrew G. Clark* *Department of Molecular Biology and Genetics, Cornell University, Ithaca, New York 14853 and †School of Biological Sciences, University of Nebraska, Lincoln, Nebraska 68588 Manuscript received February 1, 2010 Accepted for publication February 6, 2010 ABSTRACT Partial diallel crossing designs are in common use among evolutionary geneticists, as well as among plant and animal breeders. When the goal is to make statements about populations represented by a given set of lines, it is desirable to maximize the number of lines sampled given a set number of crosses among them. We propose an augmented round-robin design that accomplishes this. We develop a hierarchical Bayesian model to estimate quantitative genetic parameters from our scheme. For example, we show how to partition genetic effects into specific and general combining abilities, and the method provides estimates of heritability, dominance, and genetic correlations in the face of complex and unbalanced designs. We test our approach with simulated and real data. We show that although the models slightly overestimate genetic variances, main effects are assessed accurately and precisely. We also illustrate how our approach allows the construction of posterior distributions of combinations of parameters by calculating narrow-sense heritability and a genetic correlation between activities of two enzymes.

C

ROSSES among inbred lines (diallel crosses) have a long history in quantitative genetics (Sprague and Tatum 1942; Lynch and Walsh 1998). When one is interested in easily measurable properties of a particular, and limited, set of lines, the full diallel provides an efficient means of estimation of many parameters of interest to, for example, animal and plant breeders (Lynch and Walsh 1998). In evolutionary genetics, however, a common aim is to estimate the properties of populations the inbred lines represent. This requires a large sample of lines from a given population, and the crosses have to be replicated to control for environmental variance. Furthermore, the phenotypes assayed can be complicated or expensive to determine. Because the number of crosses in the full diallel grows very rapidly with the number of lines, this design is not practical for such purposes (Kempthorne and Curnow 1961; Lynch and Walsh 1998). A variety of designs where only a subset of possible crosses is performed have been developed over the years (Lynch and Walsh 1998). Ideally, one would like a scheme where the number of crosses grows linearly with the number of lines sampled, while maintaining the ability to estimate important parameters relevant to evolutionary processes, such as narrow-sense heritability. The round-robin design (Wayne et al. 2004), where each Supporting information is available online at http://www.genetics.org/ cgi/content/full/genetics.110.115055/DC1. 1 Corresponding author: 221 Biotechnology Bldg., Cornell University, Ithaca, NY 14853. E-mail: [email protected] Genetics 185: 361–373 (May 2010)

line participates in one cross as a male and one as female, comes close. This design is related to the circulant crosses proposed by Kempthorne and Curnow (1961). The round-robin scheme is attractive because the number of crosses is equal to the number of lines, and thus the number of lines sampled with a given number of crosses is maximized. However, this design provides low power to estimate the main effects of lines [general combining ability (Sprague and Tatum 1942; Lynch and Walsh 1998)], and for X-linked traits assayed in males the task is impossible (Wayne et al. 2004). Another possibility is to use only inbred lines themselves (e.g., Ayroles et al. 2009a). While useful for association mapping, this approach precludes separate estimates of general and specific (i.e., pertaining to particular cross) combining ability. Furthermore, the effect of inbreeding on the genetic architectures of traits is largely unknown (Lynch and Walsh 1998; Charlesworth and Charlesworth 1999). Recent studies show that mRNA levels of many genes, especially those involved in metabolism and stress response, are appreciably altered in inbred lines of Drosophila melanogaster (Kristensen et al. 2006; Ayroles et al. 2009b), potentially changing the associations of nucleotide variants with phenotypes. Inbreeding may thus confound the results of quantitative genetic experiments. While partial diallel designs are an attractive option when the number of crosses is limited, estimation of quantitative genetic parameters from such schemes is in general difficult (Lynch and Walsh 1998). Traditional

362

A. J. Greenberg et al.

ANOVA and maximum-likelihood approaches can deal with relatively simple designs, especially in the absence of replication (Lynch and Walsh 1998). Development of Bayesian methods, coupled to Markov chain Monte Carlo (MCMC) sampling, greatly expanded the flexibility to estimate parameters from crossing schemes of arbitrary complexity (Sorensen and Gianola 2002). Moreover, Bayesian posterior distributions of parameters of interest take account of uncertainty in estimates of all other elements of a model. Finally, MCMC methods allow sampling from distributions of arbitrary combinations of parameters. This feature greatly expands the set of questions that can be probed statistically. Animal breeders have long embraced Bayesian statistics (Gianola and Fernando 1986; Blasco 2001; Sorensen and Gianola 2002; Thompson et al. 2005) and these methods are popular in population genetics and quantitative trait mapping (Beaumont and Rannala 2004). In contrast, applications in evolutionary quantitative genetics are still rare (Walsh 2001; O’Hara et al. 2008). We therefore set out to explore the utility of a Bayesian approach in estimating parameters of interest to evolutionary quantitative geneticists. In particular, we examined hierarchical modeling, a method that is popular in the social sciences (Gelman and Hill 2007) but has not been widely adopted by geneticists. Our purpose is not to demonstrate the superiority of the Bayesian perspective—we believe that the choice between Bayesian and maximum-likelihood approaches should be dictated by pragmatic considerations. Rather, we seek to demonstrate that Bayesian modeling can provide an accurate and useful alternative to traditional methods, even in the face of a complicated experimental design. We propose an augmented round-robin crossing design and a hierarchical Bayesian model to estimate quantitative genetic parameters from it. This approach is motivated by our interest in assessing within- and between-population variation in a variety of metabolic traits, building on the work of Clark (1989). We augment the round-robin design with crosses between populations and measurements of the inbred lines lines themselves. Thus, each line participates in four crosses and is also ‘‘selfed,’’ enabling us to separate general and specific combining ability. We replicate the crosses to measure environmental effects. The nesting of environmental replicates within crosses, and lines within populations, is naturally modeled using a hierarchical Bayesian approach (Gelman et al. 2004; Gelman and Hill 2007). We implement such a scheme and perform simulations to show that it performs well in estimating quantitative genetic parameters of interest to us. We illustrate the utility of our approach with an example using real data. The techniques we describe should be useful in a variety of settings in evolutionary quantitative genetics, as well as in animal and plant breeding.

METHODS

Simulated data: The crossing scheme we simulated is depicted in Figure 1A. We deterministically set mean values for four populations (population means: A ¼ 8.0, B ¼ 9.0, C ¼ 8.5, and D ¼ 6.0) and generated 15 line samples from a normal distribution with the mean equal to the corresponding population mean. We ran two sets of simulations—with high and low narrow-sense heritability (h2) (Lynch and Walsh 1998). We list the parameters of these simulations in Table 1. The standard deviations (SD) for the normal distributions we used to generate line means were 1.5 for the high- and 0.5 for the low-heritability sets. To simulate crosses, we drew values from normal distributions with means equal to midparental values and standard deviations 1.0 (high heritability) and 0.4 (low heritability). Cross values were then adjusted for inbreeding and outcrossing by subtracting a value drawn from a normal distribution with mean 2.0 and SD 0.25 for inbred lines and adding a value drawn from a normal distribution with mean 1.0 and SD 0.1 for the between-population crosses. In keeping with the real data sets we are generating in the lab, we hierarchically simulated two kinds of environmental effects (Figure 1B). First, we generated two ‘‘blocks’’ of replicates by drawing two values for each cross from a normal distribution with mean equal to the cross mean and SD 0.7 (high heritability) and 1.0 (low heritability). Next, we generated ‘‘replicates’’ from blocks by drawing three values for each block mean from a normal distribution with SD 0.5 (high heritability) and 1.5 (low heritability). In the absence of epistasis, narrow-sense heritability can be estimated by the formula h2 ¼

ðs2A

s2A 1 s2D 1 s2E Þ

(Lynch and Walsh 1998). For crosses of inbred lines, s2A can be estimated by calculating s2GCA ðs2A ¼ 2s2GCA Þ and s2D by s2SCA ðs2D ¼ s2SCA Þ, assuming epistatic interactions are weak (Lynch and Walsh 1998). s2A (additive genetic variance) is the variance among line means from a given population. In our case, the true s2A is the square of the SD of the distribution from which our simulation drew the line means for each population (see supporting information, File S1 for details). Dominance is the deviation of each cross from the mean of the two parental lines. Thus, the true s2D in our simulations is the square of the SD of the distribution used to pick cross values from parental line means. The environmental variance, s2E , is the sum of the block and replicate variances. Plugging in the values for these parameters enumerated above, h2 for high- and lowheritability cases is 0.56 and 0.07, respectively. In our experiments, we measure enzyme kinetics on each replicate. We do this by providing a substrate to fly

Bayesian Quantitative-Genetic Model

extracts and measuring changes in optical density of a substrate or product of an enzymatic reaction with time. The slope of the resulting line is the maximum enzyme rate, Vmax. This rate is thus measured with some error, which is in practice variable from reaction to reaction. To investigate the effect of this on the inferences of model parameters, we simulated Vmax values from normal distributions centered on the corresponding replicate values. To mimic the real world variation in assay quality, standard deviations for these distributions were drawn from truncated normal distributions (constraining SD to be above zero) with mean 20 times the replicate mean and SD 4 times the replicate value. We then generated points from these distributions and reestimated the means (replicate values) and their standard deviations. In measuring enzyme activity, some reactions run aberrantly and produce outliers. To assess the impact of outlier observations on inference, we duplicated each data set, randomly drew 1% of observations, and multiplied each by a value drawn from a uniform distribution between 3 and 6. The regression standard deviations were multiplied by the same number. Thus, each data set has a ‘‘regular’’ version and a version with outliers. Each simulated data set consists of 1080 sets of estimated enzyme rates with their standard deviations. We implemented the simulations in R (R Development Core Team 2008) and generated 500 data sets for each of the heritability values. The R script we used, with detailed annotations, is included in File S1. Real data: To check the performance of our model on real data, we chose two well-characterized enzymes from the pentose phosphate pathway, glucose-6phosphate dehydrogenase (G6PD, EC 1.1.1.49) and 6phosphogluconate dehydrogenase (6PGD, EC 1.1.1.44). We performed the assays as described in Clark and Keith (1989). We chose 92 D. melanogaster lines from five populations and inbred them for 12 generations. Nineteen lines came from The Netherlands [courtesy of Z. Bochdanovits (Bochdanovits and Jong 2003)], Ithaca, New York (collected in 2004 by E. M. HillBurns and B. P. Lazzaro), and Tasmania (sent in 2003 courtesy of A. A. Hoffman); 17 were from Beijing [provided by C. Aquadro (Begun and Aquadro 1995)]; and 18 were from Zimbabwe [ZH lines are from Harare and ZS lines are from Sengawa, provided by C. Aquadro (Begun and Aquadro 1993), and ZW lines are from Victoria Falls, provided by W. Ballard]. The names of the lines are listed in the data file (enz_data.tsv) in File S4. We crossed the lines in a scheme similar to the one depicted in Figure 1A and replicated the crosses in two blocks, with three replicates for each block. Five male flies were collected per replicate, weighed, crushed, and resuspended in buffer (Clark and Keith 1989). We performed the kinetic assays in 96well plates, each plate containing extracts from one replicate of a particular kind of cross (i.e., inbred lines

363

or within- or between-population crosses). The complete data set is provided in the enz_data.tsv file in File S4. Normal model for simulated data: We describe the important features of our Bayesian hierarchical model here. A full description of the distributions of parameters and the Gibbs sampling scheme can be found in FileS2. Traditionally, data similar to ours are analyzed using linear mixed models. The effects of, say, line are either ‘‘random,’’ if one is interested in the population the lines come from, or ‘‘fixed,’’ if one is interested in the particular sample of lines (Box and Tiao 1973), although different interpretations of these terms are also in use (Gelman 2005). Bayesians consider all parameters as random variables; thus the terms fixed and random effects are confusing. We use the terms ‘‘sample’’ and ‘‘population’’ parameters, when referring to groups of variables, or specifically identify them as means or variances. Depending on the biological question we are considering, we may be interested in either or both kinds of variables, and the strength of the Bayesian approach is that it allows us to simultaneously estimate them all (Box and Tiao 1973; Gelman 2005). Bayesian analysis of mixed-effects models, together with MCMC sampling from posterior distributions, is generally successful (Box and Tiao 1973; Sorensen and Gianola 2002; Gelman et al. 2004). However, a number of computational difficulties exist (Gilks and Roberts 1996; Gelman et al. 2004), largely due to posterior correlations between the population and sample parameters (Gelfand et al. 1995; Gilks and Roberts 1996; Gelman et al. 2004). Bayesian hierarchical models make use of the data structure to alleviate these difficulties using two basic techniques. First, the population parameters are centered on means for the corresponding level [i.e., the sample parameter values—hierarchical centering (Gelfand et al. 1995)], rather than on zero as is customary in mixed-effects models (Lynch and Walsh 1998; Sorensen and Gianola 2002). Thus, for example, in our model the block values mbl j are distributed normally with the mean among all blocks derived from the same cross (mcross k½j ) and variance s2bl . Second, for a given level in the hierarchy, the likelihood for a sample parameter depends on the information from the level below, while the prior comes from the level above (Gelman et al. 2004; Gelman and Hill 2007). To illustrate this, we consider estimation of a block mean, mbl j . The likelihood is rep

N ðmi2j ; s2rep Þ; rep

where mi2j is the mean of all the replicate values that belong to the same block j (hence the notation i 2 j), and s2rep is the variance among replicates. The prior is

364

A. J. Greenberg et al. 2 N ðmcross k½j ; sbl Þ;

mbl j N

where mcross is the value for cross k that the block bek½j longs to, and s2bl is the among-block variance. Given that the normal prior is conjugate to the normal likelihood, the posterior distribution of mbl j is rep

mbl j

N

ðnj =s2rep Þmi2j 1 ð1=s2bl Þmcross k½j nj =s2rep 1 ð1=s2bl Þ

ln mln $l½j 1 m#l½j

1 ; nj =s2rep 1 1=s2bl

!

ð1Þ (Box and Tiao 1973; Gelman et al. 2004; Gelman and Hill 2007), where nj is the number of replicates in each block j and the other parameters are as before. Thus, each value mbl j is pulled toward the relevant cross mean, and the strength of the pulling is dependent on the precision of the estimate of mbl j from the data on one hand and the precision of the estimate of the cross mean on the other (Gelman et al. 2004; Gelman and Hill 2007). Note also that unbalanced designs (where nj are not the same across blocks, for example) are easily dealt with in this framework. The grand mean does not have a level above it, and therefore we assign it an improper flat prior in the model. It has been shown that with the types of model we are considering here, this choice of prior for the grand mean still leads to a proper posterior (Gelman et al. 2004). An additional complication is that this parameter is estimated from a small number of populations and thus the assumption of normality seems unrealistic. Instead, we assumed that population means come from a Student’s t distribution with 3 d.f. To estimate line means, we consider all the crosses a given line participates in. Because there are three types of crosses—the selfed inbred lines and the within- and between-population crosses—we first have to correct for cross type effects. In principle, this is straightforward. Taking the within-population crosses as the baseline, we regress the cross means on two indicator variables: one takes the value of 1 when the cross is a selfed inbred line and 0 otherwise; the other takes the value of 1 only when the cross is between populations. The intercept is then the specific combining ability, or the cross mean after correcting for cross type. Since the coefficients of inbreeding and outcrossing (binbr and boutc) are not further modeled (and are given flat priors), this is the well-known fixed-slope, variable-intercept regression (Gelman and Hill 2007). However, in practice it turns out that this approach leads to posterior correlations between the regression coefficients and the amongblock variance (Gelman et al. 2004), leading to poor mixing and often to convergence of s2bl estimates to zero (not shown). To alleviate this problem, we implemented parameter expansion (Gelman et al. 2004, 2008) when estimating the cross means. To do this, we modeled block means as

2

! SCA 1 X b 1 bk½j ; s2bl ;

ln where mln $l ½j and m#l ½j are the line means for the female and male parents of the cross, X is the matrix of indexes for inbreeding and outcrossing, b ¼ ðbinbr ; boutc ÞT , and SCA is the coefficient for specific combining ability, bk½j which we model as normally distributed with mean 0 and ln SCA variance s2SCA . Note that ðmln $l ½j 1 m#l ½j Þ=2 1 X b 1 bk½j ¼ cross mk½j of Equation 1. Thus, we abandon hierarchical centering for this level and use a mixed model. To break SCA 2 the posterior correlation among b, bk½j , and sbl , we SCA introduce parameter expansion by replacing bk½j with a product of two values:

SCA ¼ aSCA gk½j bk½j

gk N ð0; s2g Þ: The three values, aSCA, gk, and sg are not independently identified in the model and thus do not converge on any value (Gelman et al. 2004, 2008). However, the parameters we are interested in can be recovered at each iteration of the Markov chain: ¼ mSCA k ¼

ln mln $l½k 1 m#l½k

2 ln m$l½k 1 mln #l½k 2

1 bkSCA 1 aSCA gk

and s2SCA ¼ ðaSCA Þ2 s2g : These variables are defined in the model and converge on the appropriate values, but because they are products of undefined parameters that are free to wander in the sample space, the posterior correlations are broken, leading to faster and more accurate convergence (Gelman et al. 2004, 2008). In our implementation of the Markov chain, we update the b-vector as a batch by first calculating the mean of mbl values that belong to the same cross k, ln subtracting the current values for ðmln $l ½k 1 m#l ½k Þ=2 and SCA a gk, and regressing the resulting value on X without an intercept as suggested in section 18.4 of Gelman and Hill (2007). Once the inbreeding and outcrossing effects have been controlled for, we use the values of mSCA to calculate line effects or general combining abilities. For each line l, we examine all the crosses it participates in (five in total). We take the mSCA k2l value of the selfed cross as is. For the other crosses, we subtract the line mean for the opposite parent from twice the mSCA k2l value for each cross. Thus, we have five estimates of the line mean from five different crosses. The likelihood for the line effect is a normal distribution with mean equal to average of the


five estimates and variance s2SCA . The variance among line effects is twice the general combining ability variance (2 3 s2GCA ). Likelihoods for estimates of variances are inverse-x2. We allow the three kinds of crosses to have different s2rep and s2bl , on the basis of a priori biological considerations and some preliminary analyses of real data. We assume a single value for both s2SCA and s2GCA . We used flat improper priors on s [s() } 1 (Gelman 2006; Gelman and Hill 2007)]. We also tried proper Gamma priors, but these tend to shrink variance estimates to zero (Gelman 2006; Gelman and Hill 2007), leading to slow mixing and often to convergence failure (not shown). Each replicate value in our data set is the slope of the corresponding kinetic curve or the enzyme maximal rate (Vmax). We treated these two ways. First, we simply took the point estimates of Vmax. Second, we modeled Vmax values hierarchically with the prior being the mean value for the relevant block of replicates and the likelihood based on estimates of replicate values and their standard deviations computed as described at the end of the Simulated data section. Student’s t models for outliers: For data sets with outliers, constructed as discussed above, we implemented a further extension of our model. Instead of a normal, we used Student’s t distributions for replicate means: rep

mi

2 tnrep ðmbl j½i ; srep Þ:

Student’s t distributions have fatter tails than the normal, and therefore estimates of mbl j should be more robust to outliers (Gelman et al. 2004). For a given value of the variance, samples generated from a t distribution are expected to contain more observations that are far from the mean compared to samples generated from a normal. Thus, outlier observations do not inflate estimates of the variance, and the extent to which this is true depends on the degrees of freedom. As the degrees of freedom increase, a t distribution approaches the normal with the corresponding mean and variance. We report analyses of the simulated data sets using t distributions with 3 and 6 d.f. Three degrees of freedom is the smallest value for which the mean and the variance for the t distribution are finite. We implemented sampling from t distributions on the basis of their interpretation as mixtures of normals (Gelman et al. 2004). To improve mixing, we used parameter expansion for variance estimates (as suggested in section 11.8 of Gelman et al. 2004). Details of the sampling scheme are in File S2. Computation: Because we assume normal or t distributions for sample parameters throughout, the posterior distributions for all variables are available in closed form. Thus, we can take advantage of the efficiency of Gibbs sampling (Gilks et al. 1996) to construct the marginal posterior distributions for all parameters. We used

365

blockwise updating for variables of the same level, with the exception of line effects, which had to be updated one at a time. The full sampling schemes for both the normal and the Student’s t model can be found in File S2. To set initial values, we first calculated approximate point estimates of parameters. For example, block means are means of the corresponding replicates, cross means are means of blocks, and so on. We then picked starting values from overdispersed distributions centered on these approximate estimates. For example, starting values for block means were picked from normal distributions with means equal to the approximate block estimates and variances 4 times the approximate estimate of s2rep . We picked initial values for standard errors from uniform distributions with lower bounds 0.2 times the approximate estimate and upper bounds 5 times the estimate. We implemented the sampling algorithms in R. We provide an example R script that implements the Student’s t model, with detailed annotations, in File S2. We ran three independent chains for each data set, with 500 iterations of burn-in followed by 1000 sampling iterations. We processed the chains using the coda package in R (version 0.13-4) (Plummer et al. 2009). If we suspected lack of convergence, we looked at timeseries graphs of parameters. In a number of cases, lack of convergence was clearly due to insufficient burn-in because an initial value happened to be far from the eventual estimate. In these cases, we reran all three chains with new initial values. For each variable, we compiled several statistics across the 500 high-heritability and 500 low-heritability data sets. For each data set, we noted whether the true value fell within the posterior 95% credibility interval, the fractional difference (i.e., the absolute difference divided by the true value) between the median of the posterior distribution and the true value, the time-series estimate of the coefficient of variation (CVR) (the standard error divided by the true value; time-series SE calculated by the summary function in coda) of the posterior distribution, and the Gelman–Rubin convergence diagnostic (Gelman and Rubin 1992, as implemented in coda). Restricted maximum-likelihood estimates of parameters: As an alternative to our Bayesian hierarchical models, we used restricted maximum likelihood (REML) to estimate parameters using a traditional mixed-effects model. We employed the REML implementation provided by the lmer() function in the lme4 package in R (Bates and Maechler 2009). Gelman and Hill (2007) provide several examples of the use of this function for the analysis of hierarchical models. We used point estimates of Vmax as the response variable and treated the replicate, block, cross (SCA), and line (GCA) effects as random. Population and cross-type effects were fixed. We describe the details of our model, the R code, and the results in File S3.

366


Model and computation for real data: To analyze real kinetic data for G6PD and 6PGD, we implemented a version of the Student’s t model with 3 d.f. described above. The model was programmed in R. We provide the annotated script, along with the raw data, in File S4. Three complications arise in the real data that we did not incorporate into our simulations. First, a small subset (0.4%) of the data is missing. Second, because the number of lines in each population is unequal, the design is somewhat unbalanced. Third, we performed the enzyme assays in 96-well plates and inspection of slope estimates clearly indicated the presence of plate effects. In the analyses presented here, we did not implement missing data imputation, treating the missing data as not collected. Thus, the problem reduces to slight imbalance in the replication structure. The model detailed above does not assume equal sample sizes at any level of the hierarchy, and thus we did not have to make any changes in our algorithms to accommodate the missing data and the imbalance in the crossing design. We corrected for plate effects within the Gibbs sampler, as described in File S3. We ran the Gibbs sampler for both enzymes simultaneously, updating variables at each level in blocks as described above for simulated data, one enzyme at a time. We generated five chains, sampling 20,000 iterations per chain after a 2000-iteration burn-in. We then thinned the output, retaining 2000 values per chain. RESULTS

Cross design and analysis: Our choice of cross design and analysis methods is motivated by our empirical research. We study evolution of metabolic function within and between populations, using D. melanogaster as a model organism. Looking at within-population variation, we are interested in the potential for adaptation to new conditions (and hence in narrow-sense heritability), as well as the effects of deleterious mutations that are maintained by mutation–selection balance (and hence in the effects of inbreeding). Between-population variation should be informative about local adaptation, and thus we want to compare population means and look for potential genetic incompatibility by examining between-population crosses. The regulatory network that binds enzymes together can be probed by measuring genetic correlations among enzyme activities (Clark 1989). For this, we need estimates of line effects (general combining abilities). Given these considerations, we developed a modified round-robin crossing design (Figure 1A). It combines round-robin (Wayne et al. 2004) crosses within populations, where each line participates in one cross as a male and in another cross as a female. To this, we added the inbred lines themselves and crosses between

Figure 1.—Data structure. (A) The diallel table. Solid and shaded squares mark the crosses that we performed. The table is subdivided by population. The order of the lines is from left to right for the female axis and from top to bottom for the male axis; i.e., the top left corner of the table represents the cross of the A1 male to the A1 female. Within each population, the F1 round-robin crosses between lines (off-diagonal) are depicted as solid squares, whereas the inbred lines themselves (diagonal) are shown as shaded squares. (B) Hierarchical levels of the data. For the ‘‘crosses’’ level, the solid arrow represents a ‘‘selfed’’ inbred line, dashed arrows represent within-population crosses, and dotted arrows show betweenpopulation crosses.

populations. Each line participates in two interpopulation crosses, once as a male and once as a female, each with a different population. This makes the betweenpopulation crosses balanced both at the line level (all lines participate in an equal number of crosses) and at the population level (each population is crossed to all other populations an equal number of times). Like the regular round-robin crossing design, the number of crosses grows linearly with the number of lines evaluated, allowing adequate sampling to estimate population parameters. Since each line is crossed five times, we can additionally separate specific and general combining ability. To estimate environmental effects, we perform structured replications of our crosses (Figure 1B). Each cross is carried out on two batches of food (block level) and three samples are taken from each batch (replicate level). This allows us to capture two kinds of random environmental effects: variations in food quality from batch to batch and fluctuations within vials where we rear the flies. We then perform enzyme kinetic assays on each replicate, obtaining estimates of enzyme maximal rates with some error (Clark 1989; Clark and Keith 1989). As is clear from Figure 1B, our data are structured hierarchically. Each level is completely nested in the level above, with the exception of the cross level that is not perfectly nested within lines (two different lines participate in within- and between-population crosses) or populations. To make use of this structure in the data, it is natural to employ hierarchical models (Gelman and Hill 2007) to estimate the parameters of interest. Furthermore, a Bayesian approach is attractive because it allows for models of arbitrary complexity and permits


a full account of uncertainty in all parameters (Sorensen and Gianola 2002; Gelman et al. 2004; Gelman and Hill 2007). However, the limitation of this approach is that the posterior distributions of variables are not available, and thus one has to resort to numerical methods to estimate them (MCMC in this case). While in theory MCMC methods should accurately estimate posterior distributions given an infinite chain (Gilks et al. 1996), in practice one has to demonstrate that they do so after a reasonable number of iterations. Therefore, we generated simulated data and used them to test the performance of our models. Analysis of simulated data: We constructed simulated data sets to answer three main questions. First, we wanted to know if our combination of crossing design and modeling approach can be used to estimate parameters of interest to us and to quantify the accuracy and precision of such estimates. Second, we wished to know whether taking full account of uncertainty in enzyme activity estimates was necessary or if we could save on computational time and use point estimates of slope values. Third, because high-throughput processing of enzymatic assays yields a small but significant number of outliers, we were interested in gauging their effect on inference and testing strategies to minimize the influence of outliers. To answer these questions, we generated two sets of simulated crosses with environmental replication: one with low narrow-sense heritability (h2 ¼ 0.07) and one with high (h2 ¼ 0.56; see methods for details). We also added outlier observations to each group of data sets. We chose values of fixed parameters (Table 1) in the simulations to fall in the range we typically see with real data. We generated 500 data sets with each type of simulation and analyzed them under three models. The first used point estimates of kinetic curve slopes and assumed normal distributions of replicates. The second also assumed normal distributions for replicates, but modeled the uncertainty in their estimation. For data sets with outliers, we also used a third model that assumed Student’s t distributions for replicates. Because t distributions have fatter tails than the normal distribution, we expect estimates of sample parameters from this model to be robust to outliers (Gelman et al. 2004). On the basis of a priori biological considerations, we modeled environmental variances (s2rep and s2bl ) separately for each type of cross, although they were the same in the simulations. The results we show are for the variances from the within-population crosses. The results for the variances from other cross types were indistinguishable. We first turn to estimates of population parameters. For technical reasons, it is easier to estimate standard deviations rather than variances when programming in R, and this has an added advantage of bringing the values for these parameters to the same scale as sample variables. Population parameters are notoriously diffi-

367 TABLE 1

Parameters common to all simulations Low h2

High h2

Population parameters A B C D

8.0 9.0 8.5 6.0

8.0 9.0 8.5 6.0 Cross type parameters

b outc b

inbr

2.0 1.0

srep sbl sd sa

1.5 1.0 0.4 0.5

2.0 1.0 Standard deviations 0.5 0.7 1.0 1.5

cult to estimate accurately (Lynch and Walsh 1998; Sorensen and Gianola 2002; O’Hara et al. 2008). We quantified the quality of our estimates several ways. To determine the accuracy of our models, we plotted the range of differences between the estimated medians of posterior distributions and the true values (Figure 2, top panels). Precision of the estimates is reflected in the extent of the spread of their posterior distributions, quantified by standard errors of the sampling Markov chains (Figure 2, bottom panels). We scaled these statistics by true values to enable comparisons among variables. Another measure of accuracy is the fraction of the time the true value falls within the posterior 95% credibility interval (95% C.I.) (Table 2). Formally, under the Bayesian paradigm this statistic does not carry the same straightforward interpretation as for frequentist confidence intervals (Box and Tiao 1973). However, from a practical standpoint this measure still reflects the accuracy of parameter estimation, as well as the degree to which the uncertainty of the estimate is captured by a model. Finally, we assessed convergence properties of the various models using the Gelman–Rubin diagnostic (Gelman and Rubin 1992). Values of this statistic close to 1.0 indicate that the particular set of chains has converged. We considered a particular run to have failed to converge if the value of the diagnostic was .1.5 and counted the number of such instances for every model–data set combination (listed as percentage of failure to converge in Table 2). We tried a range of cutoff values (not shown). Our conclusions are not sensitive to the particular value chosen. Overall, our model and cross combination yields estimates of population parameters that are close to reality (Figures 2 and 3). When outliers are absent, full treatment of uncertainty in slope estimates does not make an appreciable difference, with the exception of among-replicate standard deviations (Figure 2). In the latter case, however, the model yields overestimates that

368


Figure 2.—Accuracy and precision of standard deviation estimates. For each parameter, we show a pair of graphs. The top one plots the fractional difference between medians of estimated posterior distributions of a given parameter and true values (that is, the difference divided by the true value). The bottom plot shows coefficients of variation of the estimated posterior distributions (time-series SE divided by the true value; see methods). The box plots represent data across 500 simulated data sets with high heritability. Each plot represents results under six scenarios: 1 and 2 are for simulations without outliers and 3–6 are for those with outliers (see methods). 1 and 3 were analyzed using point estimates of slopes; 2 and 4 were analyzed modeling the uncertainty in slope values and assuming that replicates are normally distributed; and 5 and 6 were analyzed modeling the slopes, but with t3 (5) and t6 (6) distributions for replicates.

are falsely precise. In contrast, in the presence of outliers the model with point estimates yields results that are both incorrect and imprecise (Figure 2 and left set of columns in Table 2). Standard deviations are generally overestimated (with the exception of sbl), leading to drastic underestimation of heritability when true h2 is high (model 3 in Figure 3A). Furthermore, the model fails to converge more frequently than any other data set–model combination (Table 2). The deficiencies of the model with point estimates can be corrected to a large degree by modeling the uncertainty of slope estimates (model 4 in Figure 2 and Table 2). Most of the estimates (with the exception of sbl) are further improved by using t models for replicates (models 5 and 6 in Figure 2). However, in every case our models overestimate the genetic standard deviations(sSCA and sgca) by 20%. This does not result in misestimates of h2 when true heritability is high (Figure 3A), although for data sets with low heritability, where these patterns are similar but exaggerated, we tend to overestimate h2 by almost twofold (Figure 3B). Modeling the slopes also improves the probability of convergence (right set of columns in Table 2), although problems still persist for estimates of sbl when heritability is low. Convergence problems are essentially eliminated with Student’s t models (columns 5 and 6 in the right set of Table 2).

Despite problems with estimating population parameters, our models do well when evaluating sample variables (Figure 4 for high-heritability data sets; the results for low-heritability data sets and for population means in both groups of data sets are indistinguishable). For data sets with outliers, the model with point estimates still produces unbiased estimates, but they are significantly less precise. However, although the credibility intervals are wider (model 3 in Figure 4), they are only slightly more likely to include the true observation (Table 2). We can improve the precision of estimates by modeling the slopes. Even more precise estimates can be obtained with t models. As is the case with standard deviations, the model with point estimates often fails to converge when outliers are present. Modeling slopes alleviates the problem, and using t models results in further improvements (Table 2). As we detail in the methods section, posterior distributions of sample parameters at a given level in the hierarchy incorporate information from the higher levels through prior distributions. Thus, if our models perform correctly, extreme observations at a given level should be pulled toward the values for corresponding variables above them. For example, line effects of lines with true means below the corresponding population means should be overestimated. Conversely, line means for lines with values higher than their population values should be underestimated. Furthermore, the


369

TABLE 2 Credibility interval coverage and convergence rates % in posterior 95% C.I. Parameter srep sbl sSCA sgca h2 binbr boutc mline srep sbl sSCA sgca h2 binbr boutc mline

% convergence failure

1

2

3

4

5

0.0 96.0 86.6 70.8 96.2 93.0 96.4 88.6

81.8 94.8 87.4 71.6 85.6 93.2 97.0 88.7

0.0 97.0 33.0 47.6 11.2 96.2 97.2 89.2

18.8 95.0 81.8 67.8 78.6 93.8 96.6 88.5

23.8 6.4 90.4 63.4 93.8 93.2 95.4 89.6 0.2 6.2 89.6 0.0 1.4 95.6 96.0 92.2

72.0 95.8 28.0 0.8 5.8 96.4 98.8 91.5

95.6 95.6 36.0 1.0 2.2 96.2 98.8 91.3

2.0 96.6 0.4 0.4 65.2 98.0 98.4 92.9

28.6 94.6 12.0 0.4 22.2 96.6 97.2 92.5

6

1

2

3

4

5

6

h2 ¼ 0.56 60.0 0.0 12.0 0.0 88.0 0.2 62.4 0.0 89.8 0.0 94.2 0.0 94.8 0.4 89.1 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.4 36.6 13.0 3.2 2.8 24.6 24.2 14.2

0.0 6.2 3.2 0.2 0.2 4.0 5.8 2.0

0.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.8 0.4 0.4 0.4 0.2 0.0 0.2 0.2

h2 ¼ 0.07 28.6 0.4 18.4 1.8 79.6 0.6 0.0 0.6 2.4 0.6 95.2 0.6 95.6 0.8 91.8 0.5

0.4 2.6 1.8 0.0 0.2 1.0 2.0 0.6

1.0 24.8 9.4 4.2 3.6 12.4 12.0 6.7

0.0 11.0 3.2 1.2 0.6 3.4 3.4 1.2

0.0 0.0 0.2 0.2 0.2 0.0 0.2 0.2

0.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0

Models 1–6 are as in Figure 2. The set of columns on the left gives fractions of all simulations where the true value of a given parameter lies in the posterior 95% credibility interval of its estimate. The right set of columns reflects how often the Gelman– Rubin diagnostic (see methods) exceeds 1.5. Other threshold values lead to the same conclusions. Data for line means of all lines were pooled.

more extreme the true value is for a line, the more severe the discrepancy, with opposite signs for values that are low or high. To see if this indeed occurs, we plotted differences between estimated and true values of line effects against differences between true line effects and corresponding population means (botttom right panel of Figure 4; both differences scaled by the corresponding true population mean). As expected, we see an appreciable negative correlation between the two variables (Pearson’s r ¼ 0.18). It is not our main purpose to compare maximumlikelihood and Bayesian methods of inference. Nevertheless, we analyzed our simulated data sets using REML estimates of a traditional mixed-effects model (see File S3 and methods for details). As expected, we find that

for estimates of environmental standard deviations and fixed effects, this model behaves similarly to our normal model with point estimates of slopes. However, the genetic effect variances (specific and general combining ability) are frequently underestimated, especially when heritability is low (see the figures in File S3), leading to underestimates of heritability. This problem is particularly severe when heritability is low and outliers are present. Furthermore, when heritability is low, the estimates of sgca, and particularly sSCA, become unstable, with deviations from true values highly dependent on the particular data set. Analysis of real kinetic data: Simulated data sets are very useful when assessing model accuracy because the correct results are known a priori. However, simulations

Figure 3.—Accuracy and precision of narrow-sense heritability estimates. The plots are arranged and labeled as in Figure 2. Estimates are from highheritability (A) and lowheritability (B) simulations.

370


Figure 4.—Accuracy and precision of estimates of sample parameters. The box plots are arranged as in Figure 2. For line means (mline), we pooled information for all 60 lines across 500 simulated data sets. The scatterplot shows the relationship of fractional differences between true line means (mline) and the corresponding true population means (mpop) on the x-axis and the fractional difference between the posterior estimates of line means (m ˆ line ) and the corresponding true line means (mline) on the y-axis.

are necessarily idealized, and although we did our best to make them realistic (by, for example, including outliers), it is impossible to simulate every conceivable scenario. We therefore wanted to assess the performance of our model on a real data set. We analyzed data for two well-characterized enzymes in the pentose phosphate pathway, G6PD and 6PGD of D. melanogaster. Assays for these enzymes are well established (Clark and Keith 1989) and a number of studies have documented genetic variation in their activities (e.g., Bijlsma 1980; Wilton et al. 1982; Clark 1989). We used a version of our Student’s t model with 3 d.f. (model 5 discussed above), but with correction for assay–plate effects (see methods for details). The choice of model was dictated by the presence of outliers in data sets for both enzymes. Furthermore, the outliers inflated the estimates of srep to the extent that problems arose with evaluation of block standard deviations: estimates of sbl came close to zero, and this created numerical errors and lack of mixing in the posterior distributions of other parameters. We present estimates of selected parameters in Table 3. Because previous studies employed only inbred lines and different replication strategies to control for environmental effects, it is difficult to make direct comparisons between our results and preceding findings. However, we can say that the relative magnitudes of srep and sbl we see are broadly similar to those reported by Clark (1989), although we find a larger genetic component of variation in G6PD and 6PGD than in that

study. Furthermore, we find that the magnitudes of environmental standard deviations vary among cross types (Table 3). One way to assess the quality of our line effect estimates is to calculate the genetic correlation between the activities of the two enzymes. A positive correlation has been repeatedly documented before (Bijlsma 1980; Wilton et al. 1982; Clark 1989). Since we were estimating the parameters for each enzyme simultaneously in our Gibbs sampler, we were able to sample from the distribution of the correlation between line effects. We indeed see the expected positive genetic correlation (Table 3). Strikingly, our estimate is almost identical to the one obtained by Clark (1989): his value is 0.31 (0.18, 0.43), whereas ours is 0.30 (0.17, 0.42). DISCUSSION

We set out to develop a combination of a sparse partial diallel crossing scheme and modeling approach to estimate a number of parameters of interest in evolutionary quantitative genetics. For evolutionary applications, we are interested in learning about populationlevel processes, and therefore large samples of lines are required. The full diallel, which incorporates crosses among all lines, even when excluding reciprocal crosses and inbred lines, is not ideal for such applications because the number of crosses grows quickly with the number of lines assessed (Kempthorne and Curnow 1961; Lynch and Walsh 1998). We settled on a modified


371

TABLE 3 Analysis of G6PD and 6PGD kinetic data G6PD Parameter s1rep s2rep s3rep s1bl s2bl s3bl sSCA sgca h2 binbr boutc rG6PD;6PGD b

a

Estimate 0.69 0.77 0.28 0.61 0.68 0.56 0.24 0.33 0.20 0.26 0.42

6PGD 95% C.I.

Estimate

(0.62, 0.75) 0.24 (0.70, 0.85) 0.19 (0.25, 0.32) 0.16 (0.52, 0.73) 0.29 (0.55, 0.83) 0.14 (0.49, 0.64) 0.18 (0.12, 0.33) 0.11 (0.28, 0.41) 0.12 (0.15, 0.27) 0.17 (0.09, 0.44) 0.17 (0.58, 0.26) 0.55 0.30 (0.17, 0.42)

95% C.I. (0.21, 0.28) (0.17, 0.21) (0.14, 0.18) (0.22, 0.59) (0.10, 0.17) (0.15, 0.21) (0.08, 0.15) (0.10, 0.15) (0.09, 0.24) (0.24, 0.09) (0.62, 0.47)

a

Environmental standard deviation indexes are as follows: 1, within-population crosses; 2, inbred lines; 3, between-population crosses. Narrow-sense heritability is calculated with s1rep and s1bl . b Genetic correlation (between-line effects).

round-robin design that incorporates inbred lines and crosses within and between populations. Each line is selfed and also crossed four times. Thus, it is possible to partition genetic variation into specific and general combining abilities (Sprague and Tatum 1942; Lynch and Walsh 1998) and still have the number of crosses grow linearly with the number of lines. While in principle this crossing design should be useful for our purposes, it is not straightforward to analyze using traditional methods. An important stumbling block is the inclusion of three different types of crosses, making it necessary to control for cross-type effects when estimating SCA. Other difficulties include measurement errors for phenotypes (in our case, enzyme activities), two levels of environmental replication, and slight imbalances in the crossing scheme caused by unequal sampling of lines from populations. Furthermore, we are interested in evaluating distributions of deterministic combinations of parameters (for example, narrow-sense heritability and genetic correlations between enzyme activities). Bayesian approaches can fulfill these requirements (Blasco 2001; Sorensen and Gianola 2002) and have been extensively used in animal breeding and quantitative trait mapping (Gianola and Fernando 1986; Sorensen and Gianola 2002; Beaumont and Rannala 2004). Taking advantage of the structure of our data, we implemented a set of hierarchical Bayesian models (Gelman et al. 2004; Gelman and Hill 2007). We made extensive use of recently developed computational techniques—parameter expansion (Gelman et al. 2004, 2008) and hierarchical centering (Gelfand et al. 1995; Gilks and Roberts 1996)—to speed and improve convergence. Our emphasis is on estimating parameters given a model, rather than significance testing. Although model comparison techniques (Raftery 1996)

and methods for integrating over various models (Phillips and Smith 1996) are available, in our case the structure of the data and biological considerations drive the selection of the basic modeling framework. We consider three types of models that differ only in the statistical approach to dealing with replicates and choose among these approaches on the basis of their performance on simulated data. We used four simulated data sets: two with low heritability and with and without outliers among slope values and two with high heritability (also with and without outliers). We then assessed precision and accuracy of estimates from our models, compared with true values. As expected (Gelman 2005), estimates of sample-level parameters (e.g., line effects and cross-type coefficients) are more accurate than evaluations of population-level variables (e.g., standard deviations and heritability). In particular, genetic standard deviations tend to be overestimated by our models, although narrow-sense heritability is still accurately estimated when it is high. This is consistent with an analysis of an approximate circulant cross of Scotts pine data (Waldmann and Ericsson 2006). This data set had only one level of replication, one type of cross, and one population. Thus, overestimation of genetic effects is probably not due to the complexity of our experimental scheme. We find our combination of crossing design and model can effectively partition genetic components of variation and adequately estimate and control for cross-type effects. The latter is particularly important, because large inbreeding and outcrossing effects are present in the real data sets we seek to model (Table 3). Overall, we find that when no outliers are present taking point estimates of slopes instead of modeling them does not appreciably affect inference. However, in the presence of outliers, modeling the slopes is essential.

372


Further improvements in sample parameter estimates can be achieved when we use Student’s t models for replicates. While analyzing empirical data, we further found that t models are often necessary to prevent Markov chains from getting stuck on estimates of sbl close to zero, a behavior that leads to numerical problems (values too small for computers to handle) and failure of proper mixing (Markov chains not moving adequately through the values that belong to the posterior distribution of a given parameter). As a comparison, we used REML to estimate the parameters of a mixed-effects model to analyze our simulated data sets. We found that in general, this approach gives results that are similar to our normal model with point estimates of slopes. However, in contrast to our Bayesian estimates, the genetic variances are underestimated by the mixed-effects model. This is particularly pronounced when outliers are present and the heritability is low. In this case, estimates of sgca and particularly sSCA become fairly unstable, with the deviations from true values highly dependent on the data set (see File S3). This problem does not appear to arise in our Bayesian estimates. An important feature of the Bayesian approach is the emphasis on parameter estimation over hypothesis testing. In the traditional non-Bayesian framework, one would, for example, test the significance of the additive genetic effect. If the effect is significant, the line means are estimated separately, and if it is not, line effects are pooled to estimate, say, population parameters. In the hierarchical Bayesian framework, we use the relative information from the population and line levels to partially pool line effects, with the degree of pooling determined by the data (Gelman and Hill 2007). We illustrate this behavior in the bottom right plot of Figure 4, where we see that estimates of means for lines with extreme values are shrunk to their respective population means. Another example of this approach is our separate treatment of environmental variances for each cross type. Thus, instead of testing for the significance of the gene-by-environment interaction and seeking a biological interpretation if it is significant, we are able to ask a biologically driven question while building our model. We find differences among some of the cross types (Table 3). However, even when differences are not ‘‘significant’’ (for example, posterior intervals for srep of 6PGD in within- and between-population crosses overlap substantially), we do not eliminate the extra parameters from the model. While this treatment does not capture all the possible kinds of gene-by-environment interactions, it provides an example of a biologically driven attempt to dissect them. Because it is easy to construct posterior distributions of any deterministic combination of parameters, such tests can be driven by biological interest rather than statistical convenience. We provide examples of esti-

mates of narrow-sense heritability and genetic correlations between activities of two enzymes: G6PD and 6PGD. We show that our approach recovers the wellknown positive correlation between these two enzymes (Bijlsma 1980; Wilton et al. 1982; Clark 1989) that is likely due to the toxicity of an intermediate compound (Hughes and Lucchesi 1977, 1978). Because the previous studies used only inbred lines to calculate this correlation, concerns have been raised (Zera and Harshman 2001) that this type of observation is an artifact of making recessive deleterious alleles homozygous through inbreeding. Since we estimate our correlation values from general combining abilities on the basis of heterozygous as well as inbred lines, we can exclude this caveat. The combination of crossing scheme and modeling approach we propose provides a powerful framework for estimating parameters important in evolutionary quantitative genetics. Although our primary motivation was the assessment of natural variation in metabolic attributes including enzyme activities, it can be readily extended to any situation when phenotypes are measured with error and are laborious or expensive to assess. We thank B. Logsdon and two anonymous reviewers for helpful comments on the manuscript and A. Coventry for help with computational methods. This research was funded by a National Institutes of Health grant to A.G.C. and L.H.

LITERATURE CITED Ayroles, J. F., M. A. Carbone, E. A. Stone, K. W. Jordan, R. F. Lyman et al., 2009a Systems genetics of complex traits in Drosophila melanogaster. Nat. Genet. 41: 299–307. Ayroles, J. F., K. A. Hughes, K. C. Rowe, M. M. Reedy, S. L. RodriguezZas et al., 2009b A genomewide assessment of inbreeding depression—gene number, function, and mode of action. Conserv. Biol. 23: 920–930. Bates, D., and M. Maechler, 2009 lme4: Linear Mixed-Effects Models Using S4 Classes. Beaumont, M. A., and B. Rannala, 2004 The Bayesian revolution in genetics. Nat. Rev. Genet. 5: 251–261. Begun, D. J., and C. F. Aquadro, 1993 African and North American populations of Drosophila melanogaster are very different at the DNA level. Nature 365: 548–550. Begun, D. J., and C. F. Aquadro, 1995 Molecular variation at the vermilion locus in geographically diverse populations of Drosophila melanogaster and D. simulans. Genetics 140: 1019–1032. Bijlsma, R., 1980 Polymorphism at the G6pd and 6Pgd loci in Drosophila melanogaster. IV. Genetic factors modifying enzyme activity. Biochem. Genet. 18: 699–715. Blasco, A., 2001 The Bayesian controversy in animal breeding. J. Anim. Sci. 79: 2023–2046. Bochdanovits, Z., and G. Jong, 2003 Temperature dependent larval resource allocation shaping adult body size in Drosophila melanogaster. J. Evol. Biol. 16: 1159–1167. Box, G. E. P., and G. C. Tiao, 1973 Bayesian Inference in Statistical Analysis. Wiley Classics, New York. Charlesworth, B., and D. Charlesworth, 1999 The genetic basis of inbreeding depression. Genet. Res. 74: 329–340. Clark, A. G., 1989 Causes and consequences of variation in energy storage in Drosophila melanogaster. Genetics 123: 131–144. Clark, A. G., and L. E. Keith, 1989 Rapid enzyme kinetic assays of individual Drosophila and comparisons of field-caught D. melanogaster and D. simulans. Biochem. Genet. 27: 263–277.

Bayesian Quantitative-Genetic Model Gelfand, A. E., S. K. Sahu and B. P. Carlin, 1995 Efficient parametrisations for normal linear mixed models. Biometrika 82: 479–488. Gelman, A., 2005 Analysis of variance: why it is more important than ever. Ann. Stat. 33: 1–31. Gelman, A., 2006 Prior distributions for variance parameters in hierarchical models. Bayesian Anal. 1: 514–534. Gelman, A., and J. Hill, 2007 Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, Cambridge, UK. Gelman, A., and D. B. Rubin, 1992 Inference from iterative simulation using multiple sequences. Stat. Sci. 7: 457–472. Gelman, A., J. B. Carlin, H. S. Stern and D. B. Rubin, 2004 Bayesian Data Analysis, Ed. 2. CRC Press, London. Gelman, A., D. A. van Dyk, Z. Huang and W. J. Boscardin, 2008 Using redundant parameterizations to fit hierarchical models. J. Comput. Graph. Stat. 17: 95–122. Gianola, D., and R. L. Fernando, 1986 Bayesian methods in animal breeding theory. J. Anim. Sci. 63: 217–244. Gilks, W. R., and G. O. Roberts, 1996 Strategies for improving MCMC, pp. 89–114 in Markov Chain Monte Carlo in Practice, edited by W. R. Gilks, S. Richardson and D. J. Spiegelhalter. Chapman & Hall, London. Gilks, W.R.,S.RichardsonandD.J. Spiegelhalter,1996 Introducing Markov chain Monte Carlo, pp. 1–20 in Markov Chain Monte Carlo in Practice, edited by W. R. Gilks, S. Richardson and D. J. Spiegelhalter. Chapman & Hall, London. Hughes, M. B., and J. C. Lucchesi, 1977 Genetic rescue of a lethal ‘‘null’’ activity allele of 6-phosphogluconate dehydrogenase in Drosophila melanogaster. Science 196: 1114–1115. Hughes, M. B., and J. C. Lucchesi, 1978 Dietary rescue of a lethal ‘‘null’’ activity allele of 6-phosphogluconate dehydrogenase in Drosophila melanogaster. Biochem. Genet. 16: 469–475. Kempthorne, O., and R. N. Curnow, 1961 The partial diallel cross. Biometrics 17: 229–250. Kristensen, T. N., P. Sørensen, K. S. Pedersen, M. Kruhøffer and V. Loeschcke, 2006 Inbreeding by environmental interactions affect gene expression in Drosophila melanogaster. Genetics 173: 1329–1336. Lynch, M., and B. Walsh, 1998 Genetics and Analysis of Quantitative Traits. Sinauer Associates, Sunderland, MA.

373

O’Hara, R. B., J. M. Cano, O. Ovaskainen, C. Teplitsky and J. S. Alho, 2008 Bayesian approaches in evolutionary quantitative genetics. J. Evol. Biol. 21: 949–957. Phillips, D. B., and A. F. M. Smith, 1996 Bayesian model comparison via jump diffusions, pp. 215–240 in Markov Chain Monte Carlo in Practice, edited by W. R. Gilks, S. Richardson and D. J. Spiegelhalter. Chapman & Hall, London/New York. Plummer, M., N. Best, K. Cowles and K. Vines, 2009 coda: Output Analysis and Diagnostics for MCMC. R Development Core Team, 2008 R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. Raftery, A. E., 1996 Hypothesis testing and model selection, pp. 163–188 in Markov Chain Monte Carlo in Practice, edited by W. R. Gilks, S. Richardson and D. J. Spiegelhalter. Chapman & Hall, London/New York. Sorensen, D., and D. Gianola, 2002 Likelihood, Bayesian and MCMC Methods in Quantitative Genetics. Springer, New York. Sprague, G. F., and L. A. Tatum, 1942 General vs. specific combining ability in single crosses of corn. J. Am. Soc. Agron. 34: 923–932. Thompson, R., S. Brotherstone and I. M. White, 2005 Estimation of quantitative genetic parameters. Philos. Trans. R. Soc. B 360: 1469–1477. Waldmann, P., and T. Ericsson, 2006 Comparison of REML and Gibbs sampling estimates of multi-trait genetic parameters in Scots pine. Theor. Appl. Genet. 112: 1441–1451. Walsh, B., 2001 Quantitative genetics in the age of genomics. Theor. Popul. Biol. 59: 175–184. Wayne, M. L., Y.-J. Pan, S. V. Nuzhdin and L. M. McIntyre, 2004 Additivity and trans-acting effects on gene expression in male Drosophila simulans. Genetics 168: 1413–1420. Wilton, A. N., C. C. Laurie-Ahlberg, T. H. Emigh and J. W. Curtsinger, 1982 Naturally occurring enzyme activity variation in Drosophila melanogaster. II. Relationships among enzymes. Genetics 102: 207–221. Zera, A. J., and L. G. Harshman, 2001 The physiology of life history trade-offs in animals. Annu. Rev. Ecol. Syst. 32: 95–126.

Communicating editor: K. W. Broman

Supporting Information http://www.genetics.org/cgi/content/full/genetics.110.115055/DC1

A Hierarchical Bayesian Model for a Novel Sparse Partial Diallel Crossing Design Anthony J. Greenberg, Sean R. Hackett, Lawrence G. Harshman and Andrew G. Clark

Copyright © 2010 by the Genetics Society of America DOI: 10.1534/genetics.110.115055

2 SI 

A. J Greenberg et al.

FILE S1 Generating Simulated Data File S1 is available for download as a compressed file (.zip) at http://www.genetics.org/cgi/content/full/genetics.110.115055/DC1. This file includes: simulation.Rnw simulation.pdf


3 SI 

FILE S2 Student-t Models for the Analysis of Simulated Data File S2 is available for download as a compressed file (.zip) at http://www.genetics.org/cgi/content/full/genetics.110.115055/DC1. This file includes: sim_analysis.Rnw sim_analysis.pdf models.pdf

4 SI 


REML AnalysisFILE ofS3Simulated Data REML Analysis of Simulated Data

A. J. Greenberg, S. R. Hackett, L. G. Harshman, and A. G. Clark

As a comparison with our Bayesian approach, we ran a simple traditional mixed-effects model. We define the model as follows: ♀ln pop ln bl Vmax(plkji) = µ♂ + µp♀pop + β inbr xinbr + β outc xoutc + b♂ + bpl + bsca p k k pl plk + bplkj + "plkji , where the random effects are 2 "plkji ∼ N(0, σrep ) 2 bbl plkj ∼ N(0, σbl )

2 bsca plk ∼ N(0, σsca ) ln 2 ∼ N(0, σ♂ ) b♂ pl ln ♀ln 2 bpl ∼ N(0, σ♀ln )

2 2 2 = σ♀ σa2 ≈ 2σgca ln + σ♂ln

♀pop pop The effects of the male parent’s (µ♂ ) and female parent’s (µp ) population are treated p inbr as fixed, as well as the effects of inbreeding (β ) and outcrossing (β outc ). The notation follows that for the rest of the paper. The effects of the male and female parent are assumed independent. Therefore, twice the general combining ability variance is the sum of variances 2 among male and female parental lines (σ 2 and σ♀ ln ), and is equal to additive genetic ♂ln 2 variance (σa ) if epistasis is ignored. Note also that the integers p and l index population and line pairs respectively. To estimate the parameters of the model, we used the lmer() function from the lme4 package in R (see main text for references). We load a simulated data set (generated with code from Supplement 1) and run the function, redirecting the summary to a variable: source("sim_data.R") res