| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
ANIMAL GENETICS |
* Departamento Mejora Genética Animal, INIA, Crta. de la Coruña km 7.5, 28040 Madrid, Spain;
and
Departamento Genética, Facultad de Biología, UCM, Ciudad Universitaria s/n, 28040 Madrid, Spain; and
CSIRO Livestock Industries, Queensland Bioscience Precinct, 306 Carmody Rd., St Lucia 4067, Queensland, Australia
 |
Abstract |
|---|
 Top Abstract INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONLITERATURE CITED |
|---|
Key Words: Bayesian mixed linear model • differential gene expression • false discovery rate • microarray • model selection • normalization
|
INTRODUCTION |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONLITERATURE CITED |
|---|
). Data normalization and analysis search for methods that remove experiment-wide systematic sources of bias in the detection of differentially expressed genes, while controlling both false negative and false positive inferences. Single gene analyses of microarray data are commonly used. However, this approach ignores valuable information from other genes in the estimation and adjustment of nonbiological or technical factors and of dispersion parameters. Hierarchical models that combine population and single gene information have recently been proposed (Smyth, 2004) to circumvent this problem. A global approach that uses data from all genes in a joint analysis was also proposed as a way to provide more power of detection (Kerr et al., 2000). However, assumptions on within-gene variability may be an issue in this type of analysis (Cui et al., 2005).
Mixed linear models have been used to perform global analysis of microarray data (Reverter et al., 2003a; Byrne et al., 2005). However, little attention has been paid to model selection. Several alternative models, in terms of definition and assumptions on effects, in addition to sources of heterogeneity of variances can be investigated (Kerr et al., 2002; Altman, 2005). Other researchers have stressed the need to define the appropriate model to accommodate the experimental design, and, therefore, to account for the sources of variation that will permit the appropriate choice of error term for hypothesis testing (Churchill, 2002; Rosa et al., 2005).
The objective of this paper was to compare several candidate models for data normalization and analysis to identify differentially expressed genes in an experiment comparing gene expression between 2 skeletal muscles, in a Bayesian context.
|
MATERIALS AND METHODS |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONLITERATURE CITED |
|---|
). Materials
A bovine cDNA microarray from fat and muscle cDNA libraries (Lehnert et al., 2004) was used to study differential expression of genes between 2 skeletal muscles in cattle. The microarray consisted of 9,600 cattle probes, comprising complementary DNA sequences (cDNA) from the bovine skeletal muscle and subcutaneous fat cDNA libraries, as well as Meat Animal Research Center (Clay Center, NE) bovine expressed sequence tag libraries and oligonucleotide probes. A total of 3,411 elements had functional annotation and the remaining elements were anonymous cDNA clones from non-normalized muscle and fat-derived cDNA libraries. The main classes of genes represented on the microarray were contractile proteins, muscle metabolic enzymes, extracellular matrix genes, and lipogenic enzymes, all of which were associated with the principal cellular components of muscle tissue. The microarray also included a collection of candidate genes from the adipogenesis and protein turnover literature (Lehnert et al., 2006, 2007).
Ten male calves of the Avileña-Negra Ibérica breed were used in this experiment. Avileña-Negra Ibérica is a beef cattle breed reared under extensive conditions and in which meat quality is emphasized as a breeding objective. Calves were fattened with the same diet and had similar ages at slaughter. Tissue samples from 2 skeletal muscles, psoas major (PM) and flexor digitorum (FD), were collected at slaughter. Psoas major has been characterized as a mixed-glycolytic muscle, whereas FD was identified as oxidative (Moreno-Sánchez, 2007). Tissues were transferred into liquid nitrogen and stored at –80°C until RNA isolation. A total of 20 microarray slides were used following a loop design as shown in Figure 1. The design layout aimed at controlling dye bias while allowing for an efficient between-muscles contrast; direct comparisons were made between muscles, within and between individuals, taking into account both muscle and individual differences. The dyes were swapped between muscles to account for the effect of dye bias.
|
Microarray slides were scanned with the GenePix 4000B scanner (Axon Instruments, Union City, CA) at a resolution of 5 µm (Photo Multiplier Tube) values ranging from 550 to 700 and laser power 100%). GenePix Pro 5.0 software (Axon Instruments) was used first to obtain the Cy3 and Cy5 image [green (G) and red (R) fluorescence, respectively] TIFF files, locating the spots in the microarray with the appropriate grid, and then to acquire the intensity values from both channels at each spot.
To avoid bias in the identification of differentially expressed genes, criteria for data acquisition dealing with the relationship between signal and background and the accuracy of the spot intensity (Dudoit et al., 2000; Tran et al., 2002) were applied. Data acquisition criteria aim at reducing the noisy nature of the experimental conditions being surveyed. A FORTRAN 90 computer program was developed to visually assess the quality of the raw data and to filter invalid spots. Raw data quality was first assessed by performing MA plots (Dudoit et al., 2000). Such plots showed that all values were homogeneously distributed around zero, indicating no dependencies of log2 (R/G) or M value (differences among dyes and muscle) on the overall intensity level (A). Invalid spots were excluded for further normalization and analysis according to the following criteria: First, the spots flagged by the GenePix as "not found," "absent," or "bad" were removed. Second, the spots for which the background signal intensity exceeded the foreground signal intensity were also discarded. Third, to address signal quality, spots with signal to noise ratio <1 and spots with a mean to median correlation <0.85 were rejected (Reverter et al., 2003a). The above-mentioned criteria had to be satisfied by both channels in each spot, otherwise the information from both channels was discarded. Finally, probes that did not have at least 2 spots per array and were in more than 1 array were ignored. In total, 269,712 intensity readings from both channels were deemed as valid and represented a total of 8,538 clones.
These intensity records were background-corrected by subtracting the background from the foreground intensities in both colors and log2-transformed so that they approximated a normal distribution.
Statistical Analysis
Data normalization and analysis were performed jointly for all genes. The complete generic model is described as follows in matrix notation:
where y is a vector of measured intensities. Vector β includes effects related to technical matters that may affect intensity readings (e.g., print-tip effect, dye effect, etc.), and lacking, therefore, biological meaning. The spot location was modeled either by an effect of the array (a, with 19 levels) or the combination of the specific location of the spot within the array (array_block combination, ab, with 912 levels). Additionally, main effects included muscle (m, with 2 levels) or dye (d, with 2 levels). First-order interactions between array dye (ad, 38 levels), or array_block x dye (abd with 1,824 levels) and muscle x dye (md with 4 levels), were also considered in the most complete model. Vector g contained the clone or gene effects that described the existing variability among clones included in the analysis (8,538 levels). Spot-to-spot variability was accounted for by either vector ag (the vector of gene x array interaction with 60,998 levels) or abg (the vector of gene x array_block interaction with 73,523 levels). Vector dg is the vector of dye x gene interaction (17,076 levels). It models the clone-specific dye effects that occur when some groups of genes show more affinity for a specific fluorescent channel, regardless of the treatment. Vector mg is the muscle x gene interaction (17,076 levels), which is the biological result we aimed for. This vector contains the differences among genes within muscle, or (which is the same) the differences in gene expression between muscles. Vector cg (59,549 levels) is a vector containing the differences in gene expression across animals. This effect reflects the differences in genotype among animals, which may also explain differences in the level of expression among genes. The matrices X, Zg, Zag, Zdg, Zmg, and Zcg are the corresponding incidence matrices associated with the effects described previously. From here to the end of this section, ag refers in general to both treatments of the effect (ag or abg).
For the Bayesian analyses, data were assumed to be generated from a multivariate normal distribution, according to the following stochastic process:
where 
g2, ag2, dg2, mg2, and cg2 are the variances among clones or genes, genes within array, genes within dye channel, genes within muscle, and genes within individuals, respectively. Finally, Ri is the error variance matrix containing the residual variance associated with each observation. Ri refers to the within-gene variance. The subindex i indicates that residual variances may change across the levels of i. Several assumptions were made relative to this component. Assumptions are described in the following section.
Conjugate priors were used for location and dispersion parameters. Proper-vague priors were employed to assign a low degree of belief to the prior information and to allow the computation of the marginal density of the data and the Bayes factor (BF). The BF is not well defined for improper priors (Gelman et al., 1995). For the location parameters, N distributions were assumed as follows:
where 
is a positive, scalar hyper parameter, large enough to approximate a flat distribution (more precisely, = 106). Under these specifications, effects in β approach a uniform distribution, which would result in effects named as fixed in a non-Bayesian context. The remaining effects can be regarded as random effects in the non-Bayesian statistical tradition. Effects with nearly flat and normal distributions will be called fixed and random effects, respectively, herein, to allow for a shorter and more widely used notation in the microarray literature.
For the dispersion parameters, scaled inverse chi-squared (
–2) distributions were used:
where parameters 
i and Si are interpreted as degrees of belief and a priori values for the variances of the g, ag, dg, mg, cg, and e effects. A value of 4 was assigned to the degrees of belief for the –2 distributions to provide low weight to the prior information. Values for the scalars 
, Sag2, Sdg2, Smg2, Scg2, and 
were REML estimates obtained under the complete model. Values of scalars Sag2 varied according to the definition of the effect, ag or abg.
Posterior marginal inferences on parameters of interest were drawn from their corresponding conditional posterior distributions through a Gibbs sampling scheme. All conditional distributions had known forms as described in Sorensen and Gianola (2005). The specifications for the Gibbs sampling implementation were based on the Raftery and Lewis (1992) criterion, which was determined using the public domain program Gibbsit v.2.0 (Raftery and Lewis, 1995). Alternatively, the coupling method analysis (Johnson, 1996; García-Cortés et al., 1997) was performed. Finally, 100,000 iterations after a burn-in period of 20,000 rounds were carried out in all analyses.
Model Comparison
Because of the relatively large number of effects considered, an exhaustive stepwise procedure could not be performed. A full model (FM) containing all mentioned main effects and their first-order interactions was initially fitted. In the FM model, the gene effect and its interactions were considered as random effects, and the rest of the effects were treated as fixed, according to the previously specified prior distributions. Subsequently, several reduced models were analyzed. A first group of reduced models, SP_i (i = 1,3), dealt with spatial arrangements of spots within and across arrays. Thus, the spot effect was modeled either by an effect of the array (a) or the combination of specific location of the spot within the array (array_block combination, ab) and interactions involving gene x array (ag) or gene x array_block (abg). A second group of reduced models, LB_i (i = 1,2), dealt with interaction effects relative to labeling, such as dye x muscle (md) and array_block x dye (abd). Third, 4 additional models, GI_i (i = 1,4), were run to account for effects related to specific gene interactions. In a second set of comparisons, the best model from the first step was further investigated to compare the inclusion of the effects of array and gene as random or fixed effects. Two additional models were considered, one with both array_block and gene effects treated as random factors, ABr_Gr, and another model where array_block was assumed to be random and the gene effect was fixed, ABr_Gf. To do this, priors were changed according to the corresponding assumption. All models were compared under the assumption of homogeneity and heterogeneity of residual variances according to level of intensity. Six classes of intensities were considered with upper limits of 5, 7, 9, 11, 13, and 16 (in the log2 scale). A description of the various models can be found in Table 1. In addition to the heterogeneous variances linked to intensity level, treatment (muscle, in our case) and a combination of treatment and intensity level were additionally investigated as sources of heteroscedasticity for LB_1, GI_1, and ABr_Gf models. Upper limits of intensities were 5, 7, 9, 12, and 16 for each muscle.
|
; Kass and Raftery, 1995) and the cross-validation predictive densities of the data (Gelfand et al., 1992). The BF for 2 competing models, Mi and Mj, was computed as the ratio of their 2 corresponding marginal densities (MD) of the data under each model. The MD for model k was obtained from the Newton and Raftery (1994) estimator:
where H is the total number of iterations after the burn-in period, and f(y |
k(i)) is the conditional distribution of the data given the parameters involved in model k at iteration i. This estimator converges to the correct value of f(y|Mk) as H tends to infinity (Kass and Raftery, 1995).
The cross-validation predictive densities are defined as the set of univariate densities (Gelfand et al., 1992; Gelfand, 1996), f(Yr|y(r)), with r = 1, ..., n, where y(r) is the vector of observations with observation yr excluded, and Yr is a random variable with distribution f(Yr|y(r)), that gives the values of the unobserved Yr that would be more likely to be observed when the model has been fitted to y(r) (Gelfand, 1996). Model comparison was done by computing the expected value of the checking function, g = Yr – yr, for all of the data with respect to its univariate predictive density,
the best model being the one that has minimum
An importance sampling, using the joint posterior density of the parameters given the data, p(|y), as an importance distribution, was implemented to evaluate [1]. In this manner, dr can be obtained for all of the data directly from the outputs of the Gibbs sampling implemented to draw the marginal posterior densities of the parameters (Gelfand et al., 1992). The log of the MD (logMD) and D statistics [2] were obtained for each case within the Gibbs sampling process (details of the implementation can be found in López-Romero et al., 2003).
Detection of Differentially Expressed Genes
From the selected model in the previous step, both under the homogeneous and heterogeneous residual variances, estimates of the differences in expression for each gene (dg) were then obtained from the posterior means corresponding to the muscle x gene interaction (
g, PM, g, FD):
Using this definition of difference, it is assumed that estimates of all levels of gene x muscle interaction are obtained with the same accuracy.
To identify differentially expressed genes, a model-based cluster analysis was applied to the differences in gene expression between both muscles (dg). Under this approach, it is assumed that the data to be clustered come from several components, each of them with a distinct normal distribution. That is, each data point in our vector of dg is taken to be a realization from a normal mixture distribution of k (finite but unknown) components with the probability density function:
where 
(dg; µj;Vj) denotes the normal density function with mean µj and (co)variance matrix Vj and the mixing proportions 
j are constrained to be non-negative and sum to unity.
This analysis aimed at grouping clones according to dg and providing to each data point a posterior probability of being assigned to each of the significant normal components. The computer package BAYESMIX (Reverter et al., 2003b) was used. The program provides posterior probabilities of belonging to each of the components for each clone and statistical criteria (Akaike information criterion and Bayesian information criterion) that allow determination of the number of components or clusters that yield the best fit to the data. A maximum of 5 components was tested. Total number of iterations was 12,000 and the burn-in period consisted of the first 2,000 rounds. Clones were assigned to the cluster with the largest posterior probability. Estimates of false discovery rate (FDR) associated with the results coming from the best models were obtained from the mixture models as in McLachlan et al. (2005) transforming the n finite numbers of clusters into a 2-component mixture, one containing genes that are identified as differentially expressed and another group formed by the nondifferentially expressed genes. Thus,
where 
circ;0(wj) is the posterior probability that the jth gene from the subgroup of differentially expressed genes (Nr), is not differentially expressed. This has been termed local FDR by Efron and Tibshirani (2002) and may be seen as an empirical Bayes version of the Benjamini and Hochberg (1995) methodology (Efron, 2004).
|
RESULTS AND DISCUSSION |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONLITERATURE CITED |
|---|
provides a model description together with model comparison criteria, logMD, to compute BF and D statistics for each analyzed model with and without heterogeneity of residual variances due to level of intensity. The lowest values of |logMD| and D indicate the best models. Table 2 shows estimates of variance components under the assumption of heteroscedasticity. The estimates of residual variances correspond to a weighted average of the estimates in each interval according to the number of observations.
|
The FM, SP_1, SP_2, and SP_3 models were compared to study the spot effect. This effect comes from differences between the print-tips on the slide, as well as spot-to-spot variability explained by its interaction with gene effects. Including the ab effect resulted in more favorable values of both model comparison criteria compared with including the a effect (FM vs. SP_2 and SP_1 vs. SP_3). Moreover, substituting the abg with the ag interaction (SP_1 vs. SP_2 and FM vs. SP_3) had a negative and much greater effect than the ab by a substitution on both model comparison criteria. The large contribution of the abg effect was also observed for model GI_1, in which this effect was excluded. Overall, models excluding the abg interaction or substituting this factor with the ag interaction showed the worst values for both comparison criteria. Regardless of the assumption on residual variances (homogeneity vs. heterogeneity), the same pattern was observed. Spot-to-spot variation has been highlighted by several researchers as one of the major sources of bias when detecting differentially expressed genes (Wolfinger et al., 2001; Kerr et al., 2002; Byrne et al., 2005). However, in those cases, spot-to-spot variability was only related to spots across arrays and not across array_block combinations. From the results obtained in our study, it seems clear that, where we have slides with technical replication (designed to have at least 2 spots per probe) developed under noncommercial conditions, the block within array provides a significantly better descriptor of the position of the spot and should be considered instead of array to remove spot-to-spot variation.
To quantify the differences among models, 2*log(BF) was calculated (not shown). Although the BF accounts for the number of variables and the amount of data available, allowing the selection of more parsimonious models (Sorensen and Gianola, 2005), models including effects with a very large number of levels were always superior according to this criterion. However, the predictive ability statistic (D) was always consistent with the BF for the set of models considered in this step. All models compared until now performed better under the assumption of heterogeneity of residual variances. Therefore, the model of choice according to the model comparison statistics at this point was the heteroscedastic FM model. This model had the best goodness of fit without compromising its predictive ability.
The estimated variance components under the SP_i models, shown in Table 2, corroborated the results from the model comparison statistics. The amount of variation explained by the spatial arrangement of clones was not negligible, regardless of the model. The proportion of estimates of the variance associated with the interaction between the gene and the spatial location (a or ab effect) to the total estimated variation ranged from 5% for model SP_1 and model SP_3, which included the ag effect, to 12% for models FM and SP_2, which included the abg effect. On average, estimates of residual variances were affected accordingly. Thus, the estimated residual variance was one-third of that associated to model SP_1 or SP_3. This result confirms that an array_block x gene interaction better accommodates the existing spot-to-spot variability in our data, in agreement with what the criteria for model selection had been indicating.
To detect differentially expressed genes, other important sources of bias are related to dye effects mainly because of labeling efficiencies. Two different interaction effects were checked in our models. These effects were muscle x dye and array_block x dye. Results are shown in Table 1. When the md interaction was removed from model LB_1, the corresponding values of absolute value of logMD and D decreased as compared with model FM. Therefore, the interaction of dye by muscle seems to be a negligible source of bias in our data. These results also indicate that gene x muscle interaction effects that account for differences among genes within muscles (the variable of interest) are free of the bias because of the unbalance between red and green intensities across muscles. As in the previous case, 2*log(BF) was calculated, and we found that LB_1 was superior to model FM. It is important to note that, at this stage, both model comparison criteria were also consistent in the results regardless of the assumptions relative to residual variances. Furthermore, comparing models LB_2 and FM allows us to make inferences about the dye x array_block interaction (abd). Both logMD and D indicated that the FM model was better than LB_2, which was the model in which array_block x dye interaction (abd) was ignored. As a consequence, array_block x dye interaction was further maintained in our model selection procedure. Its effect pointed out the importance of the disequilibrium between red and green labeling across arrays.
The appropriateness of including interactions involving the gene and the other main factors was ascertained by the model comparison statistics in Table 1 and the estimated variance components in Table 2. The best model according to logMD and D was again LB_1, which included all gene interactions and was the most complete model at this stage. Removing either array_block x gene interaction (model GI_1) or gene x muscle interaction (model GI_3) had the greatest negative impact on the model selection criteria, as well as on the estimates of residual variances. The estimated variances for the gene x array_block and the gene x muscle interactions accounted for 12 and 4% of the total variation, respectively. Consequently, residual variances were overestimated when those interactions were excluded from the model. On the other hand, the existence of an important source of variation due to the gene x muscle interaction indicates that differential expression between muscles exists. Effects related to the RNA present in the slide are represented by the gene effect in addition to its interactions. Although the estimate of the component of muscle x gene interaction was quite consistent across models, the inclusion of nuisance effects may have an effect on the inference about differentially expressed genes through the effect on the accuracy of the estimation of differences. The effect of the model fitted on inferences about differentially expressed genes is expected to derive from differences in the accuracy of the estimation of the mg effect. Thus, overfitted models are expected to provide estimates with less accuracy and shrink toward the mean. This is because a variable number of parameters have to be estimated from the same amount of information.
In another group of comparisons, assumptions relative to the fixed or random condition of the effects of array_block and of the gene effect were compared following the literature on microarray data analyses. Thus, in the linear model context, where a mixed model analysis is applied, the array effect is often accommodated as a random effect and gene and gene x treatment interaction as fixed effects (Kerr et al., 2000; Wolfinger et al., 2001; Pfister-Genskow et al., 2005). Under this approach, it is acknowledged that spot-to-spot variation is due to a stochastic process and is not a predefined source of variation. On the other hand, other researchers (Reverter et al., 2003a; Byrne et al., 2005) consider the array effect as fixed and the gene or clone effect as random. Under this assumption, genes present in a microarray are conceptually considered as a random sample of all possible genes explaining the biological differences (Altman, 2005). As we have pointed out, in the Bayesian analysis context, this is equivalent to changing the prior information relative to location and scale parameters. We explored 3 possible combinations, the one presented in FM or LB_1, in which array_block is a fixed effect and the gene effect and its interactions are random, another model where array_block is a random effect and the gene and gene interactions with dye and muscle are fixed, and the last model where both effects and the interactions are random. According to logMD and D in Table 1, the best model was the one that accommodated the array_block as a systematic effect, and the gene or clone effect was included as a random effect (LB_1). The 2log*(BF) indicated that model LB_1 was better than the model that is most often used in the microarray literature (ABr_Gf) and was first introduced in an ANOVA context by Kerr et al. (2002). As previously found, both selection criteria identified LB_1 as the model of choice. Estimates of variance components under those models are also shown in Table 2. The total variation is reduced or increased as a consequence of model assumptions. The amount of variation attributed to the array_block effect was not negligible in any case, varying on average from 7 to 28% of the total variation in models ABr_Gf and ABr_Gr, respectively. The magnitude of the effect is close to that found in other experiments using noncommercial cDNA platforms (Pfister-Genskow et al., 2005). However, such variation was an extra source of variation, increasing the total estimated variance (from 3.73 to 4.01 for LB_1 and ABr_Gr models, respectively). The corresponding estimates for the other effects were consistent with those obtained under the reference model, which in this case was LB_1. The decision of whether to treat some effects as fixed or random may have a substantial effect on the results. According to Kerr (2003), the optimum method would be to remove systematic effects through the normalization part of the model and then include random effects. In a Bayesian context, the comparison among these types of models is related to the comparison of appropriate priors relative to the amount of information available for estimation. For the so-called fixed effects, a nearly flat prior is used, which is less informative than that assigned under the random effect situation. The best model considered array_block as a systematic effect and gene and its interactions as random effects. These results point out the differential amount of information available to estimate the effects and favor the use of informative priors, which provide more certainty to the effect with less information, the gene effect in our analysis. In summary, the LB_1 model taking into account the heterogeneity of variance according to level of intensity showed the best goodness of fit (logMD) at this stage without major changes in the predictive ability of the model (D) when compared with the homogenous residual model.
Fitting heterogeneous residual variance caused a reduction in the magnitude of the estimated variance components, except for the gene x muscle interaction (results not shown) in all models. The observed reduction varied according to the model and the component, ranging from 3 to 10% for all components. Consequently, the variance of the gene x muscle interaction relative to the total variance improved from 3 to 8% for LB_1 and GI_1 models, respectively. Such gain is expected to be at the cost of favoring the finding of genes as differentially expressed when in reality they are not for the model with the worst goodness of fit and predictive ability (GI_1). Cui and Churchill (2003) highlighted the risk of generating false positive or negative differentially expressed genes when common variance assumption is not true.
In addition to models accommodating heterogeneous variance according to level of intensity, other sources of heterogeneity were investigated. Results for model comparison and variance component estimates can be found in Tables 1 and 2, respectively. Results are shown for LB_1 with heterogeneity assumed between treatments (LB_1M) and a combination of treatment_level of intensity (LB_1MI). None of the assumptions improved the goodness of fit (logMD) or the predictive ability (D) of models relative to our model of choice, LB_1 (Table 1). The estimates of variance components obtained with the model accommodating heterogeneity due to level of intensity in each muscle (LB_1MI) were similar to the ones obtained in our previous results (LB_1). Similar results were found in all tested models.
Table 3 shows the estimates of residual variance or within-gene variance in each of the scenarios studied. Heterogeneity clearly affects the low intensity classes. This is observed regardless of the model. The magnitude of the estimated residual variances for the 2 first intervals was large compared with the other groups of intensities. Similar results have been found by Dudoit et al. (2000) and Kerr et al. (2002). In addition, it can be observed that the effect that has largely affected the magnitude of the within-gene variance differences is the array_block x gene interaction. Models that ignored this effect (GI_1) or those that included the array by gene interaction only (SP_1 and SP_3) showed the greatest magnitude of estimated residual variance in the first 2 intervals of intensity and a relatively small change in the other classes. Logarithmic transformation inflated the variance when intensities were low, although it stabilized the variance at greater intensities (Lin et al., 2008), but the impact of such transformation seems to be model dependent according to our results. The ratio of estimates of variance components changed according to the interval of intensity, and, therefore, the corresponding shrinkage factor.
|
Model-based cluster analysis was used to make inference on differentially expressed genes between muscles. The measure of differences of the estimates of gene or clone effect in each of the muscles (dg) was provided by LB_1 models with and without heteroscedastic residual variances. In both cases, the fitting of mixtures revealed the model with 3 clusters as the one that best fits the data (results not shown). For the model with heterogeneous variances, the components were as follows:
For the model with homogeneous variances, the components were as follows:
Under both models, the proposed density parameters and mixing proportions were very similar. Variances tended to be larger for the case of homogeneity. The differentially expressed clones were those whose probability of being in a cluster except the one containing nondifferentially expressed genes was the largest. Thus, differentially expressed clones were 236 and 168 for LB_1 under heterogeneity and homogeneity according to level of intensity, respectively. Regardless of the case, most of the clones had the greatest probability of being assigned to the cluster of nondifferentially expressed genes. The associated estimates of FDR were similar, being 0.238 and 0.243 for the heterogeneity and homogeneity cases, respectively. These results indicate that the heterogeneity models allow us to distinguish a larger number of genes differentially expressed between muscles. The rank correlation of the estimated differences of expression among these genes was 0.97. However, when comparing the estimates of FDR at the same level of detection, that is, for the same number of differentially expressed clones, the heterogeneity model rendered a decreased FDR than the homogeneity case (0.159 and 0.243, respectively). Estimated FDR corresponding to the differentially expressed cut-off point provided by the criterion of posterior probability greater than 0.5 was considered large. Finally, the cut-off point was established for a low value of FDR, 0.05. Under this criterion the final number of clones identified in each case varied following the same pattern of what we have previously observed. Thus, in the model with heterogeneity of variances, 82 clones were identified as differentially expressed, whereas for the homogeneity case the number was 61. Again, the model that included heterogeneity of variances associated with the level of intensities allows identification of a larger number of clones. However, the number of genes (26 vs. 23) and biological functions (5 vs. 5) associated with the differentially expressed clones were similar. When we relaxed the FDR, 198 versus 166 clones were declared as differentially expressed, referring to 37 and 35 genes, respectively. Heterogeneity of variances allows better identification of the existing redundancy of the microarray, which has been previously reported for this platform (Byrne et al., 2005). Among other biological functions, differential expression was found for genes coding for metabolic enzymes, both glycolytic and oxidative. As expected, a group of 9 genes coding for glycolytic enzymes was overexpressed in PM, the most glycolytic muscle. A group of 2 genes coding for mitochondrial enzymes involved in oxidative phosphorylation pathways was also overexpressed in PM, suggesting that PM also has a marked oxidative ability together with the glycolytic one. More details can be found in Moreno-Sánchez (2007). The 2 genes that differed between the homogeneity and heterogeneity cases were AK1 (Adenylate kinase 1) and PYGM (Phophorylase, glycogen, muscle); PYGM was 1 of the 9 genes coding for glycolytic enzymes, which were overexpressed in PM. To validate our experiment, real-time reverse transcription-PCR was performed on 4 genes including AK1. All genes were found differentially expressed in the corresponding muscle. The AK1 gene was finally confirmed as differentially expressed, constituting a biological validation of our results.
The use of shrinkage estimators, such as the one used in this analysis, is expected to prevent a large FDR (Ishwaran and Rao, 2003). This is so because the estimated expression level under the treatment of interest (the muscle in our case) is shrunk toward the mean (0 in our case) when the amount of information for a particular clone is small.
Although the use of microarray data is normally an exploratory tool, unreliable results may cause invalidation of further research when ignoring genes that are falsely declared as nondifferentially expressed. According to the results obtained in our study, more adequate modeling of residual terms could help in reducing the number of false positives. This result had been previously proven by simulation by Zhang et al. (2005). Even if our results are only valid for this specific experiment, they highlight the need to be cautious when approaching data normalization and analysis to identify differentially expressed genes. Global analysis with heterogeneity of residual variances allows combining information of similar intensity levels for estimating error variances and performing adequate inferences. This is expected to produce more robust and powerful inference than modeling residuals separately for every gene (Kerr, 2003; Cui et al., 2005; Hoeschele and Li, 2005). From a computational point of view, it is more feasible to perform global linear normalization and significance tests in 2 steps (see Wolfinger et al., 2001). This approach will provide equivalent results to the 1-step procedure, when there is not strong imbalance in the data and there is a large proportion of nondifferentially expressed genes (Hoeschele and Li, 2005). On the other hand, empirical results indicate that a gene by gene analysis increases the probability of finding false positives when differential expression is assessed (Efron et al., 2001), because estimates of error variances are obtained with less information. To avoid such problems, several authors (Smyth, 2004; Cui et al., 2005; Lo and Gottardo, 2007) have proposed procedures in terms of single gene analysis borrowing information from the variances among genes to make inferences about each gene individually. In our study, large estimates of residual variances are related to low intensity levels. Residual variance is an important part of the observed variability, which was largely reduced by adjusting factors such as the gene interactions, mainly gene x array_block interaction. Accounting for other sources of heterogeneity, such as the muscle or muscle x intensity level interaction, did not seem to improve the detection of heterogeneity among residual variation. These results and the fact that a substantial proportion of clones had a relatively small number of replicates could indicate that the within-gene heterogeneity might not differ greatly from a gene by gene approach relying on population information, with population being defined and stratified by the intensity levels.
The large differences found for the model comparison criteria across models highlight the importance of defining the factors that more accurately account for experiment-wise variability to ensure correct inference on differential expression of genes. Spatial arrangements of spots, spot-to-spot variation, labeling efficiencies, and heterogeneity of variances relative to intensity have been identified as relevant sources of variation. Our results also illustrate the importance of the experimental setup to account for possible sources of bias in the detection of differentially expressed genes. Use of proper shrinkage factors that accommodate the heterogeneity of residual variance can help control the false discovery rate.
|
Footnotes |
|---|
2 Corresponding author: cdiaz{at}inia.es
Received for publication November 7, 2007. Accepted for publication October 2, 2008.
|
LITERATURE CITED |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONLITERATURE CITED |
|---|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
