Estimation and Testing of Gene Expression Heterosis

Heterosis, also known as the hybrid vigor, occurs when the mean phenotype of hybrid offspring is superior to that of its two inbred parents. The heterosis phenomenon is extensively utilized in agriculture though the molecular basis is still unknown. In an effort to understand phenotypic heterosis at the molecular level, researchers have begun to compare expression levels of thousands of genes between parental inbred lines and their hybrid offspring to search for evidence of gene expression heterosis. Standard statistical approaches for separately analyzing expression data for each gene can produce biased and highly variable estimates and unreliable tests of heterosis. To address these shortcomings, we develop a hierarchical model to borrow information across genes. Using our modeling framework, we derive empirical Bayes estimators and an inference strategy to identify gene expression heterosis. Simulation results show that our proposed method outperforms the more traditional strategy used to detect gene expression heterosis. This article has supplementary material online.


INTRODUCTION
Heterosis, or hybrid vigor, refers to the enhanced phenotype of hybrid progeny relative to their inbred parents. Taking maize as an example, the offspring from crossing the inbred lines B73 and Mo17 are taller, mature faster, and produce greater yields than their parental lines (Hallauer and Miranda 1981). Since heterosis was scientifically documented by Darwin (1876), it has been successfully manipulated to improve many species for food, feed, and fuel industries, such as rice (Yu et al. 1997), alfalfa (Riday and Brummer 2002), tomatoes (Krieger, Lippman, and Zamir 2010), and fish (Wohlfarth 1993). Despite the intensive study and successful utilization of heterosis, the basic genomic mechanisms remain unclear (Coors and Pandey 1999;Lippman and Zamir 2007). Researchers speculate that gene expression heterosis could be among the mechanisms responsible for the phenotypic heterosis (Swanson-Wagner et al. 2006;Springer and Stupar 2007).
Due to advancements in high-throughput genomics technology (such as microarray and next-generation sequencing of RNA), it is now possible to simultaneously measure and compare expression levels of thousands of genes in parental lines and their hybrid offspring to search for evidence of gene expression heterosis. It is of particular interest to test if a gene exhibits any of the following three forms of gene expression heterosis: high-parent heterosis (HPH), low-parent heterosis (LPH), or mid-parent heterosis (MPH). A gene is said to exhibit HPH if the mean expression level of the offspring is greater than the maximum of the two parental means, LPH if the mean expression level of the offspring is smaller than the minimum of the two parental means, and MPH if the mean expression level of the offspring is not equal to the average of parental means. Let i index the genotypes of the two parents (i = 1, 2) and the offspring (i = 3). Let j (j = 1, . . . , J ) index the genes, where J denotes the total number of genes under study. We use μ ij to denote the mean expression level of gene j of genotype i. Let h j = μ 3j − max{μ 1j , μ 2j }, l j = min{μ 1j , μ 2j } − μ 3j , and m j = μ 3j − (μ 1j + μ 2j )/2. With these notations, gene j exhibits HPH, LPH, or MPH if and only if h j > 0, l j > 0, or m j = 0, respectively.
Past work on estimating gene expression heterosis using microarray data (Swanson-Wagner et al. 2006;Wang et al. 2006;Bassene et al. 2010) has used separate estimates for each gene obtained by replacing population means (μ ij , i = 1, 2, 3, j = 1, . . . , J ) with corresponding sample averages. These sample average estimators of h j and l j are problematic because they are biased and tend to underestimate h j and l j (see Appendix A). Though the sample average estimator of m j is unbiased, with only a few observations for each gene in a typical microarray experiment, the sample average estimators of m j , h j , and l j may each be highly variable.
Because high-throughput technologies measure expression of hundreds of thousands of genes simultaneously, we can utilize information across genes to improve estimation and testing of gene expression heterosis for each individual gene. For gene j , we define two latent variables α j = (μ 1j − μ 2j )/2 and δ j = μ 3j − (μ 1j + μ 2j )/2. Notice that all three types of gene expression heterosis can be written as functions of |α j | and δ j , that is, h j = δ j − |α j |, l j = −|α j | − δ j , and m j = δ j . Thus, modeling of |α j | and δ j helps to develop statistical inferences for all three types of gene expression heterosis. We model α j , the half parental difference, as a draw from a mixture of a point-mass-at-0 distribution and a normal distribution. This implies that |α j | is equal to 0 with some probability π α and equal to the absolute value of a draw from a normal distribution with probability 1 − π α . The point-mass distribution in the mixture model represents the case where the parental gene expression levels are equal, whereas the normal component corresponds to genes whose expression levels differ between the two parental lines. Similarly, we model δ j , the difference between the offspring mean and the average of the parental means, with another mixture model that has normal and point-mass-at-0 component distributions. We estimate the parameters for these mixture distributions based on observed data from all genes. Under an empirical Bayes framework, we derive posterior distributions of α j and δ j and draw inferences about gene expression heterosis from estimates of these posteriors.
We compare the empirical Bayes method with the sample average method through simulation studies where datasets were generated based on real heterosis microarray experiments or hypothetical probability models. Simulation studies show that the empirical Bayes estimators of h j , l j , and m j have smaller mean square errors (MSEs) than the sample average estimators that have been used previously. Furthermore, the empirical Bayes estimators of h j and l j are less biased than the sample average estimators, and the inferences we draw using our empirical Bayes approach are superior to traditional approaches for detecting all forms of heterosis.
The remainder of the paper proceeds as follows. Section 2 presents the proposed hierarchical model in full detail. Section 3 derives the empirical Bayes estimators and inference strategy based on the framework constructed in Section 2. Section 4 summarizes analysis results of two real experiments. Section 5 presents results of several simulation studies. Section 6 summarizes our work. R code and C code for the analysis of real experiments in Section 4, the simulation studies in Section 5, and the implementation of all our algorithms is available upon request.
In the previous section, we defined α j = (μ 1j − μ 2j )/2 and δ j = μ 3j − (μ 1j + μ 2j )/2. In order to share information across genes to improve estimation of gene expression heterosis, we propose the following models (2.1)-(2.3) for α j , δ j , and the error variance σ 2 j . Suppose that and that all α j , δ j , and σ 2 j are mutually independent. The scaled inverse χ 2 model for the error variances σ 2 1 , . . . , σ 2 J given in (2.3) follows Smyth (2004). The mixture model for α j in (2.1) models the cases where parental means are equal and where parental means differ, respectively. The hyperparameter π α specifies the proportion of genes that are equally expressed between two parents. Similarly, the mixture model for δ j in (2.2) describes the cases where the mean gene expression in the offspring is equal or not to the average of two parental means. When necessary, the model (2.1)-(2.3) may be modified as needed to better capture the features of a given dataset. For example, the mixture model could include more than one normal distribution component for α j or δ j . Although all subsequent derivations are for the model specified in (2.1)-(2.3), it is straightforward to modify our proposed approach to handle more complex models.
With no loss of information about expression heterosis, the data can be summarized by the sufficient statistics α j ≡ (ȳ 1j · −ȳ 2j · )/2, δ j ≡ȳ 3j · − (ȳ 1j · +ȳ 2j · )/2, and S 2 j (j = 1, . . . , J ). Clearly, α j and δ j are the natural sample average estimators of α j and δ j , respectively. Based on the normality assumption for y ij k , the conditional distributions of α j , δ j , and S 2 j -given α j , δ j , and σ 2 j -are (2.6) By combining (2.1), (2.3), and (2.4) it follows that the marginal distribution of α j is a two-component mixture distribution, where each component density is itself an infinite mixture of normal distributions with common mean but varying variance. This marginal distribution is determined by the hyperparameters π α , μ α , σ 2 α , d 0 , and σ 2 0 . Similarly, the marginal of the distribution of δ j has an analogous form and is determined by the hyperparameters π δ , μ δ , σ 2 δ , d 0 , and σ 2 0 . Figures 1(a) and 1(b) present histograms of empirical marginal distributions and scatterplots for α j and δ j from an alfalfa experiment and a maize experiment, respectively. Each of these datasets is discussed in more detail in Section 4, but we introduce the plots here to provide some empirical support for the model described in this section. Using methods discussed in Appendix C, we obtain estimates of our model hyperparameters, and the hyperparameter estimates determine fitted marginal densities that are plotted on top of the histograms as red lines. The fitted marginal distributions adequately capture the shape of the empirical distributions. Furthermore, the lack of correlation between α j and δ j in the scatterplots supports our model assumption of independence between α j and δ j . Thus, for both datasets, the model presented in Section 2 appears to be consistent with the main features of the data illustrated in these plots.

EMPIRICAL BAYES ESTIMATION AND TESTING OF GENE EXPRESSION HETEROSIS
Obtaining estimates of our model hyperparameters is the first step in our empirical Bayes approach. We use the method of Smyth (2004) to estimate d 0 and σ 2 0 . We estimate other hyperparameters by a combined approach of the moment method and the marginal Figure 1. Scatterplots of α j vs. δ j and histograms of empirical marginal distributions of α j and δ j (j = 1, . . . , J ) based on two real heterosis experiments. The relative sizes of α j and δ j partition the two-dimensional space virtually into subsets based on the mean expression levels of two inbred parents and their hybrid offspring as shown by dashed lines. Fitted curves represent estimated marginal densities based on the assumed model described in Section 2. (a) Alfalfa dataset. B2, B5, and F1 denote the genotypes of the two parental inbred lines and the hybrid offspring, respectively. (b) Maize dataset. B73, Mo17, and F1 denote the genotypes of the two parental inbred lines and the hybrid offspring, respectively. maximum likelihood method using data from all genes. The details of our proposed approach are provided in Appendix C. Because thousands of genes in one experiment are used to obtain the estimates of the hyperparameters, we claim that adopting the usual empirical Bayes strategy (i.e., treating these unknown hyperparameters as known and equal to their estimates) does not seriously affect the performance of the inferential procedures we describe in this section. This claim is supported by simulation studies presented in Sections 4 and 5.
Once estimates of the hyperparameters have been obtained, our goal is to draw inferences regarding expression heterosis for individual genes. Based on (2.1)-(2.6), an expression for the joint posterior distribution of (α j , δ j ) given α j , δ j , and S 2 j is derived and illustrated in Appendix B. Sampling from the joint posterior distribution of (α j , δ j ) allows us to approximate the posterior distributions of h j , l j , and m j via the relationships h j = δ j − |α j |, l j = −|α j | − δ j , and m j = δ j . Based on the form of the posterior of (α j , δ j ), one common method for sampling α j and δ j is through a Markov chain Monte Carlo (MCMC) method, such as using the Metropolis-Hastings algorithm. We have developed and implemented such a Metropolis-Hastings algorithm as illustrated in the online supplement. A good approximation of the posterior distributions of h j , l j , and m j requires a large number of draws from the joint posterior distribution of (α j , δ j ) for each gene j . By using the Metropolis-Hastings algorithm, an analysis of simulated data for only 1,000 genes took around 5 hours to complete (see more details in the online supplement). Although parallelism and/or more sophisticated sampling algorithms could help to reduce the computing time, the large number of genes in a typical transcript profiling experiment motivates us to find a faster alternative.
To substantially reduce the computing requirement and maintain good approximations of the posterior distributions of h j , l j , and m j , we derive in Appendix B the approximation to the joint posterior distribution of (α j , δ j ) given by where φ(x | μ, σ 2 ) denotes the normal density with mean μ and variance σ 2 evaluated at x, (3.2e) and the probabilities P 1j , P 2j , P 3j , and P 4j sum to 1 and are defined in Appendix B. The approximation to the joint posterior distribution of α j and δ j in (3.1) is a mixture of four joint distributions, where both α j and δ j are from point-mass-at-0 as in (3.1a), δ j is from point-mass-at-0 and α j is from a normal distribution as in (3.1b), α j is from point-massat-0 and δ j is from a normal distribution as in (3.1c), and both α j and δ j are from normal distributions as in (3.1d). The approximate posterior mixture distribution combines information from prior models and empirical observations. For example, μ α j can be expressed as a weighted average of μ α (the prior mean of α j given α j = 0) and α j (an estimate of α j based on sample means), where the weight on μ α is proportional to the prior precision of α j given α j = 0 (1/σ 2 α ), and the weight on α j is proportional to an estimate of the conditional precision of α j given α j (1/ var( α j | α j )). Similarly, σ 2 α j is the inverse of the average of the precisions 1/σ 2 α and 1/ var( α j | α j ). The approximation of the joint posterior distribution in (3.1) allows us to substantially reduce the computing requirement because we no longer need to go through a large number of MCMC iterations, but can instead directly sample from either a point-mass-at-0 distribution or a normal distribution. In addition, this leads to accurate approximations of the posterior distributions of h j , l j , and m j , as demonstrated by simulation studies in Section 5 and in the online supplement.
Given the fully specified approximate posteriors of α j and δ j and plugging in estimated hyperparameters, it is straightforward to approximate posterior distributions of h j , l j , and m j by simulation. We propose to use the estimated posterior expectations h j = E(h j | α j , δ j , S 2 j ), l j = E(l j | α j , δ j , S 2 j ), and m j = E(m j | α j , δ j , S 2 j ) as point estimators for h j , l j , and m j , respectively. Tests of HPH, LPH, and MPH, respectively, for each gene j are based on the estimated posterior probabilities We also use the estimated posterior probabilities to estimate false discovery rates (FDRs) for any family of tests that involves one test per gene. The number of positives, R(c), is the number of genes declared to exhibit a type of gene expression heterosis given the cutoff c. Taking HPH as an example, , and the estimated FDR for HPH based on estimated posterior probabilities is FDR(c) = V (c)/R(c) given cutoff c. Calculations of estimated FDRs for testing LPH and MPH are similar.

ANALYSIS OF AN ALFALFA DATASET
We used our method to analyze an alfalfa dataset on gene expression in parental lines B2 and B5 and the hybrid genotype (B2×B5). The data are available in the Gene Expression Omnibus (GEO) database (Barrett et al. 2011) with series number GSE25034. Each genotype had three biological replicates measured with Affymetrix Medicago Genome Array (Platform GPL4652). The robust multiarray average (RMA) method (Irizarray et al. 2003) was used to obtain normalized expression measures for each probeset on the array. Nonalfalfa probesets associated with the bacterial genome Sinorhizobium meliloti, along with all other probesets called absent by Affymetrix microarray suite version 5 software in all samples were filtered from the dataset (McClintick and Edenberg 2006) to leave 31,865 probesets for analysis. The hyperparameters estimated from our proposed method are summarized in row 1 of Table 1. A simulation study was conducted to assess the estimation of hyperparameters. We used the estimated hyperparameter values in Table 1 as the true parameter values to simulate data for 31,865 genes based on the hierarchical model described in Sections 2 and 3. Then, we reestimated the hyperparameters using the simulated data. We repeated this procedure 1,000 times. The estimated bias and MSE in Table 1 for each hyperparameter estimator based on these 1,000 replications show that our hyperparameter estimators are reasonably accurate and precise.
For any gene j , we sample h j , l j , and m j by simulating α j and δ j from the approximate joint posterior distribution (3.1). As an example, the contour plot of 10,000 random draws of α 20 and δ 20 from the approximate joint posterior distribution of gene "AFFX-Msa-ubq11-3_at" (gene number 20) is plotted in Figure 2. This gene has been reported to be one of the polyubiquitin genes involved in directing protein recycling and related functions (Geer et al. 2010). Based on these draws, p h 20 = P (δ 20 > |α 20 | | α 20 , δ 20 , S 2 20 ) ≈ 0.998, which gives strong evidence of HPH for this gene. As described in Section 3, we can also use the estimated posterior distributions of α j and δ j to test for any given type of heterosis while controlling FDR at a specified level. For example, we color-coded points in Fig-Figure 2. Example estimated posterior distribution for a gene exhibiting significant evidence of HPH (gene "AFFX-Msa-ubq11-3_at" in the alfalfa dataset). ure 2(a) of the online supplement to highlight genes significant at approximate FDR level 0.05 when testing for HPH (red), LPH (blue), or MPH (red, blue, or green), respectively. We also used a traditional approach based on a separate analysis for each gene to analyze the alfalfa dataset. Sample average estimates and ordinary t-tests were used to identify significant evidence of heterosis. Taking HPH as an example, ifȳ 1j · ≥ȳ 2j · , then h j = y 3j · −ȳ 1j · , and the t statistic for the one-sided ordinary t-test is h j / (1/n 3 + 1/n 1 )S 2 j . Similarly, we tested for LPH using a one-sided ordinary t-test, and we tested for MPH using a two-sided ordinary t-test of m j = 0. Given the p-values from the ordinary t-tests, we controlled FDR for the sample average method using the q-value method described by Storey and Tibshirani (2003).
The numbers of genes exhibiting significant evidence of the three types of gene expression heterosis when controlling FDR at approximately 0.05 by the sample average method and the empirical Bayes method, respectively, are in Table 2. Our empirical Bayes method identifies far more significant genes than the sample average approach.

Swanson-Wagner et al. (2009) compared gene expression of maize inbred lines B73 and
Mo17 and their hybrid offspring. They studied a total of 13,999 genes in their microarray experiment with 10 biological replicates for each of the three genotypes. The dataset is downloadable in GEO with series number GSE16136.
Log-scale expression measurements were Lowess normalized within each slide and median centered. The normalized data were analyzed with our empirical Bayes method, and the estimated hyperparameters are summarized in Table 1, row 4. The simulation described in Section 4.1 was repeated for the maize results to estimate the bias and MSE of the hyperparameter estimators. The results are summarized in the last two rows of Table 1.
Based on the posterior distributions of α j and δ j , we color-coded points in Figure 2(b) of the online supplement to highlight genes significant at approximate FDR level 0.05 when testing for HPH (red), LPH (blue), or MPH (red, blue, or green), respectively. The reported numbers of genes exhibiting each of the three types of gene expression heterosis identified by the sample average method and the empirical Bayes method, respectively, are listed in Table 2 where FDR was controlled at the 0.05 level. Once again, the empirical Bayes method reported more significant genes for all three types of gene expression heterosis than the sample average method.

SIMULATION STUDY BASED ON THE ALFALFA EXPERIMENT
We simulated 100 datasets based on the hierarchical model defined by (2.1)-(2.6) using hyperparameters equal to the estimated values from the alfalfa experiment in Table 1. For each dataset, we simulated 31,865 genes (the same number of genes in the alfalfa experiment) and three biological replicates for each genotype.
We used the empirical Bayes method to estimate h j , l j , and m j for all j . For each dataset and each type of heterosis, we ranked the estimation errors from most negative to most positive, then we averaged the estimation errors of the same rank across the 100 datasets. We used the same approach for the sample average method. The box plots of averages of ranked estimation errors are plotted in Figure 3(a) for h j , Figure 3(b) for l j , and Figure 3(c) for m j . These box plots suggest that the empirical Bayes method on average has smaller ranked estimation errors than the sample average method. The box plots also show that the averages of ranked estimation errors by the empirical Bayes method have narrower interquartile ranges than the sample average method for estimating each type of heterosis. Table 3 summarizes the averaged estimation biases and MSEs across all genes in all datasets. The empirical Bayes estimators have smaller biases and MSEs than the sample average estimators for all types of heterosis. Both the plots and statistics show substantial improvement of the empirical Bayes method over the sample average method.
For each dataset, we computed the true positive rate (TPR) given a set of fixed levels of false positive rate (FPR) for testing each type of gene expression heterosis by the sample average method and the empirical Bayes method, respectively. Then, we averaged the TPRs across 100 datasets for each given level of FPR for each of the two methods. The resulting average receiver operating characteristic (ROC) curves are plotted in Figures 3(d)-3(f) for testing HPH, LPH, and MPH, respectively. We only plotted over the range of FPR between 0 and 0.05 because FPR > 0.05 is rarely of interest in practice. The ROC curves demonstrate that our proposed tests identify more true positives than the sample average method given any fixed level of FPR for testing each type of gene expression heterosis.
By the empirical Bayes method, we estimated the FDRs for testing each type of gene expression heterosis as described in Section 3. Then, for each level of estimated FDR, the true FDRs were calculated by averaging the proportions of false positives among the declared heterosis genes across 100 datasets for each type of gene expression heterosis. We plotted the estimated FDRs against the true FDRs in Figures 3(g)-3(i) for testing HPH, All results presented above and throughout the paper are based on the approximate joint posterior density in (3.1). We compared this proposed fast and approximate method with sampling from posterior distribution via the Metropolis-Hastings algorithm. Comparison results are discussed in the online supplement. In summary, we found that whereas the estimated posterior probabilities of exhibiting HPH, LPH, and MPH are very similar for both methods, our approximate method is more than 1,000 times faster than the Metropolis-Hastings approach.

SIMULATION STUDY BASED ON THE MAIZE EXPERIMENT
The estimated hyperparameters of the maize experiment were used as the true parameter values to simulate 100 microarray datasets, each with 13,999 genes (the number of genes in the maize experiment) and 10 biological replicates for each gene of each genotype.
We analyzed these 100 datasets by the empirical Bayes method and the sample average method. The estimated bias and MSE of h j , l j , and m j estimators averaged across all genes in all datasets are summarized in Table 3. Table 3 shows that the empirical Bayes estimators are more accurate and more precise than the sample average method in estimating all types of heterosis. Figure 3 of the online supplement provides box plots, ROC curves, and FDR plots for the maize simulation results that are very similar to those displayed in Figure 3 for the alfalfa simulation in Section 5.1.

SIMULATION STUDY BASED ON PROBABILITY MODELS
To further assess the performance of the proposed empirical Bayes method, we simulated data using distributions different from those proposed in (2.1) and (2.2). Specifically, we simulated α j from a mixture distribution with a point-mass-at-0 and a t distribution with a small number of degrees of freedom (2) and a noncentrality parameter (ncp) 0.01. Independently from α j , we simulated δ j from a mixture model with a point-mass-at-0 and two normal distributions N(−0.05, 0.2) and N(0, 0.2). We simulated data for 100 microarray datasets, where each dataset contains 5,000 genes with three biological replicates for each of three genotypes. Based on the estimated hyperparameters for the alfalfa experiment and the maize experiment, we set π α = 0.8, π δ = 0.6, and simulated σ 2 j from a scaled inverse χ 2 distribution with parameters d 0 = 2.8 and σ 2 0 = 0.025. Though the data were not simulated from the proposed model, our empirical Bayes estimators, compared to the sample average estimators, have substantially smaller average bias and MSE for h j and l j as shown in Table 3. Although the averaged estimated bias for m j is slightly greater than that of the sample average method, the averaged estimated MSE is reduced by the empirical Bayes method. Figure 4 of the online supplement provides box plots, ROC curves, and FDR plots (analogous to those in Figure 3 of Section 5.1) showing that the empirical Bayes method improves upon the sample average method even though the data-generating model differs from the assumptions in (2.1) and (2.2).

DISCUSSION
Gene expression heterosis is speculated to be one possible explanation for phenotypic heterosis of traits like plant height or grain yield. One natural strategy for estimation (called the sample average method in this paper) is to simply use the sample means to replace the population means when estimating the three types of gene expression heterosis. Because there are often few observations for each gene in a microarray experiment, such estimates have high standard errors. In addition, the sample average estimators for high-parent heterosis and low-parent heterosis are also biased estimators. Furthermore, the natural t-based testing strategies that accompany the sample average method yield low detection power for all forms of gene expression heterosis.
A shrinkage method based on the sample average estimators can improve inferences on gene expression heterosis by sharing information across genes. We developed hierarchical models by placing a mixture prior model on each of two latent variables. Using an empirical Bayes method, the sample average estimates of gene expression heterosis were adjusted and shrunk toward prior means estimated from the data. The extent of shrinkage was also estimated empirically based on data. Through simulation studies based on real datasets and different probability models, we demonstrated that our empirical Bayes estimators have substantially smaller bias and MSE than the sample average estimators, and the inferences for all three types of gene expression heterosis based on the posterior probabilities also yield higher TPRs given any level of FPR than the ordinary t-tests based on the sample average estimates. We also showed that using posterior probabilities of exhibiting any type of gene expression heterosis to estimate FDR yields accurate estimates of the actual FDR. Thus, the methods we have developed provide researchers with substantially improved statistical tools for studying gene expression heterosis.
The results presented in Section 4 focus on identifying individual genes that show significant evidence of expression heterosis of various types. Rather than attempting to identify individual genes, our approach can also be used to estimate global values like the proportion of all genes that exhibit a given type of heterosis. For example, the proportion of maize genes exhibiting HPH is estimated by the average posterior probability of HPH, J j =1 p h j /J = 0.122. This estimated proportion includes genes where expression in the hybrid is only slightly higher than the maximum parental expression. In some cases, scientists prefer to concentrate on large changes in expression. With our empirical Bayes approach, it is straightforward to estimate the posterior probability of h j > k for any constant k. For example, with k = log(1.5), the average posterior probability of h j > k in the maize data is 0.0006. This indicates that genes with hybrid expression (on the original scale) more than 1.5 times that of the high parent are relatively rare.
Our work has focused on the use of gene expression measurements that can be modeled, at least approximately, by linear models with normally distributed errors. This is a standard modeling approach for microarray data. Whereas there are thousands of existing microarray datasets and more generated nearly every day, next-generation sequencing of RNA (RNA-Seq) is an increasingly popular technology for obtaining gene expression measurements. At the present state of the technology, RNA-Seq data are perhaps best treated as counts and modeled with generalized linear models involving overdispersed Poisson or negative binomial distributions (see, for example, Anders and Huber 2010;Robinson, McCarthy, and Smyth 2010;Lund et al. 2012;McCarthy, Chen, and Smyth 2012). We believe that the hierarchical modeling ideas we have proposed in the linear model framework are also likely to be very useful in a generalized linear model framework for the study of gene expression heterosis using RNA-Seq data. Developing the details of such an extension is the subject of some of our ongoing and future research.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

APPENDIX A: BIAS OF THE SAMPLE AVERAGE ESTIMATORS OF HIGH AND LOW PARENT HETEROSIS
Based on the definitions in Section 2, we can rewrite the sample average estimator of h j = δ j − |α j | as δ j − | α j |. Although α j and δ j are both unbiased estimators of α j and δ j , respectively, Thus, E( h j ) = E( δ j − | α j |) < δ j − |α j | = h j . Likewise, E( l j ) = E(−| α j | − δ j ) < −|α j | − δ j = l j . Thus, the sample average estimators of h j and l j are both biased estimators that, on average, underestimate high-parent and low-parent heterosis, respectively.

Substituting (B.2) and (B.3) into (B.1) and noting that p(S
To obtain reliable statistical inferences of α j and δ j and inferences of h j , l j , and m j , we need to draw a sufficiently large sample from the posterior distribution proportional to (B.4) for each gene j . One approach is to use the Metropolis-Hastings algorithm (see the online supplement). However, due to the inefficiency of the Metropolis-Hastings algorithm and the complex structure in (B.4), obtaining a sufficiently large sample for each of the tens of thousands of genes in a typical microarray experiment requires extensive computing power. Methods, such as parallel computing, could reduce the computing time, but the total amount of required computing power remains substantial.