Abstract
Background
In DNA microarray experiments, measurements from different biological samples are often assumed to be independent and to have identical variance. For many datasets these assumptions have been shown to be invalid and typically lead to too optimistic pvalues. A method called WAME has been proposed where a variance is estimated for each sample and a covariance is estimated for each pair of samples. The current version of WAME is, however, limited to experiments with paired design, e.g. twochannel microarrays.
Results
The WAME procedure is extended to general microarray experiments, making it capable of handling both one and twochannel datasets. Two public onechannel datasets are analysed and WAME detects both unequal variances and correlations. WAME is compared to other common methods: foldchange ranking, ordinary linear model with ttests, LIMMA and weighted LIMMA. The pvalue distributions are shown to differ greatly between the examined methods. In a resamplingbased simulation study, the pvalues generated by WAME are found to be substantially more correct than the alternatives when a relatively small proportion of the genes is regulated. WAME is also shown to have higher power than the other methods. WAME is available as an Rpackage.
Conclusion
The WAME procedure is generalized and the limitation to paireddesign microarray datasets is removed. The examined other methods produce invalid pvalues in many cases, while WAME is shown to produce essentially valid pvalues when a relatively small proportion of genes is regulated. WAME is also shown to have higher power than the examined alternative methods.
Background
The DNA microarray technique involves a series of steps, from the harvesting of cells or biopsies to the preprocessing of the scanned arrays, before analysable data are obtained. During several of these steps the quality can be affected by random factors. For instance, depending on the handling of a biological sample the mRNA can be more or less degraded [1], and the celltype composition of a biopsy can be more or less representative for the tissue in question. When arrays share sources of variation the deviations from the nominal value will be correlated. For example, two arrays from sources with degraded RNA will both tend to underestimate the expression of easily degradable genes, and two biopsies with a similar and nonrepresentative celltype composition will deviate in a similar fashion from the average expression for the ideal celltype composition.
The procedure Weighted Analysis of Microarray Experiments (WAME) [2,3] introduced a model where a covariancestructure matrix common for all genes aims at catching differences in quality by differences in variances and covarying deviations by correlations between arrays. For computations of test statistics and estimators this resulted in weighting of observations according to the estimated covariancestructure matrix, giving lower weight to imprecise or positively correlated arrays.
In order for the estimation of the covariance matrix to work in the current WAME method, the measurements of most genes must only measure noise, i.e. have an expected value of zero. This is the case in experiments where pairwise logratios are observed and where few genes are differentially expressed between any of the pairwise measured conditions. In the present paper, this crucial constraint will be relaxed to only require that most genes are nondifferentially expressed between the conditions actually being compared. Thus, nonpaired experiments can be analysed, e.g. many additional ones based on onechannel microarray data. The relaxation is realised by transforming the data to remove irrelevant information in a manner yielding transformed data with expectation zero for nondifferentially expressed genes, after which the current WAME method is applied. The transformed data are shown to give equivalent tests and estimates to those of the original data, given the corresponding covariancestructure matrices.
Problem formulation and current methods
Given a microarray experiment with n arrays and m genes, we observe for each gene g an ndimensional vector X_{g }of log_{2 }transformed values measuring mRNA abundance. In WAME the vector X_{g }is assumed to have expectation μ_{g }described by a design matrix D and a genespecific parameter vector γ_{g}, typically having one dimension per studied condition. A covariancestructure matrix Σ, common for all genes, is used to model differences in quality between arrays as different variances and shared sources of variation between arrays as correlations. A genespecific variancescaling factor c_{g }is assumed to have inverse gamma prior distribution with a global shape parameter α. Conditional on c_{g }the vector X_{g }is assumed to have a normal distribution with covariance matrix c_{g}Σ. A matrix C specifies the differential expression vector δ_{g}, describing the linear combinations of the parameters that are of main interest. Formally,
and variables corresponding to different genes are assumed independent. We want to estimate the differential expression
or we want to test for differential expression
In the current version of WAME [2,3] the estimation of the covariancestructure matrix Σ is based on a temporary assumption of expectation zero, μ_{g }= 0, for all genes, which is shown to give reasonable results if the expectation is close to zero for most genes. Thus, this is a suitable assumption for data with paired observations and few regulated genes between the pairwise measured conditions.
The WAME model can be compared with the ordinary linear model (OLM) [4],
which gives rise to the ordinary t or Ftests, and with a widely used empirical Bayes model proposed in [5] and implemented in the LIMMA package [6],
The novel feature of WAME was thus the introduction of the quality modelling covariancestructure matrix Σ.
After the introduction of WAME, a weighted version of LIMMA was proposed [7], which we will refer to as wLIMMA. There, a model with arraywise variance scales but no correlations is used,
The parameters are estimated using a restricted maximumlikelihood (REML) approach.
A widely used approach is to only consider the ordinary leastsquares estimated differential expression, often referred to as the log foldchange, here abbreviated as FC, or as the average Mvalue. In the present paper, the ranking of the genes imposed by this method will be included in comparisons, when applicable.
Results
The new version of WAME
In the current version of WAME [2,3] the covariancestructure matrix Σ is estimated using a temporary assumption that μ_{g }= 0 for most genes, i.e. that the measurements of most genes consist solely of biological and technical noise. In the new version of WAME we relax this to only assume that most genes are nondifferentially expressed, i.e. δ_{g }= 0. This allows a much larger class of experimental designs and design matrices D, most notably unpaired designs.
The trick used is to transform the data and consider
where is a suitable linear estimator of μ_{g }which is unbiased under H_{0 }and which preserves the estimability of the differential expression δ_{g}, based on only the transformed data (see Methods for details). An example is (8) below where for each gene the mean value of all arrays is subtracted.
Since the transformed data contain only noise for nondifferentially expressed genes by construction, the current version of WAME can essentially be applied to the transformed data Y_{g}. As before, the covariancestructure matrix (now Σ_{Y}) and the hyperparameter α are first estimated under a provisional assumption (now δ_{g }= 0). The maximum likelihood estimates of δ_{g }and the likelihood ratio test statistics of (3) are then computed. The tests and estimators are in fact unchanged by the transformation (7), if the covariancestructure matrices for the transformed and untransformed data are known (details given in Methods). WAME is implemented as a package for the R language [8] and is available at [9].
Evaluation on real and resampled data
To investigate the properties of the new version of WAME, two real datasets are examined. Briefly, they are analysed both using WAME and the current methods described in Background. Arrayspecific weights, pvalue distributions and rankings are produced showing clear differences between the procedures, most notably in the pvalue distributions. To investigate the power of the different procedures and to look at pvalue distributions in a controlled but realistic setting, we also analyse simulated data with real noise from the studied datasets and synthetic signal.
Description of the real datasets
Two public onechannel microarray datasets are analysed. The datasets are selected from the NCBI GEO database [10] with the criteria of having unpaired design and being sufficiently large to allow for the resamplebased simulations in Resampled data below.
In the first dataset [11], biopsies were taken from the left atrium from 20 human hearts with normal sinus rythm and 10 hearts with permanent atrial fibrillation. It is here referred to as Atrium. In the second dataset [12], mechanisms in chronic obstructive pulmonary disease, COPD, were investigated by taking lung tissue biopsies from 12 smokers with mild or no emphysema and from 18 smokers with severe emphysema. In both datasets one Affymetrix HGU133A array was used for each patient. In the present paper RMA [13] is used to obtain expression measures from the raw probewise intensities. The analyses are performed using the R language and the Bioconductor framework [14].
Analysis of the real datasets
A natural parameterisation of the included datasets is to have one parameter per condition, yielding design and hypothesis matrices
Under the null hypothesis, for each gene g and array i, an unbiased estimator of the expected value of the measurement X_{ig }is obtained by the genewise mean value over all arrays from both groups. The transformation then becomes a subtraction of that mean value, cf. (7),
Note how the transformation preserves the difference in mean value between the two groups of arrays.
If the elements in X_{g }from the different arrays had in fact independent and identically distributed noise for each fixed gene g as assumed in OLM and unweighted LIMMA, the noise in Y_{g }would have equal variances for all arrays. In Figure 1 arraywise density estimates for the transformed expression values are shown. For arrays from the same condition the distributions should be identical, reflecting the combined variability of signal and noise. For unregulated genes the expectation of Y_{g }is zero, so if the assumption of few regulated genes holds the densities from all arrays should furthermore be essentially equal. Examination of Figure 1 reveals that neither of these statements are true, indicating that some variances are highly unequal.
Figure 1. Density plots. Distribution of transformed expression values, Y, for the different arrays, in the two datasets. Colourcoding according to sample variance is used for increased clarity (blue for low variance, red for high variance). Differences in variability can be noted for both datasets.
Analogously, all pairs of arrays within each condition should have a common joint distribution and when few genes are regulated all pairs of arrays should essentially have a common joint distribution with a small negative correlation of 1/(n  1). Examination of scatter plots for all pairs of arrays shows that this is clearly not the case (some obvious examples are shown in Figure 2, all pairs are included in Additional file 1 and Additional file 2).
Additional file 1. Pairwise plots of all arrays in the Atrium dataset. Transformed expression values for all arrays in the Atrium dataset. See legend of Figure 2 for details.
Format: PNG Size: 2.8MB Download file
Additional file 2. Pairwise plots of all arrays in the COPD dataset. Transformed expression values for all arrays in the COPD dataset. See legend of Figure 2 for details.
Format: PNG Size: 3.1MB Download file
Figure 2. Pairwise plots. Transformed expression values, Y_{g}, for selected pairs of arrays within the same group. Different pairs within the same group have distinctly different correlations. Upper triangle contains scatterplots. Lower triangle contains heatmaps of the corresponding twodimensional kernel density estimates, where the majority of the genes are in the red portion of the plot, revealing important trends inside the black clouds. Diagonal red clouds in the heat maps reveal correlations between arrays. Offdiagonal numbers show estimated correlations from WAME. Diagonal boxes contain sample names and weights as well as estimated variances from WAME.
As expected from the observations above, unequal variances and nonzero correlations are estimated in the analyses with WAME, giving rise to highly unequal weights in the estimates of the differential expressions (shown in Table 1 and Table 2). In fact, the sign of the weight for some arrays even get switched compared to the sign of the weight of the other arrays from the same condition. This is an effect of strong correlations combined with unequal variances. It is an issue which is further addressed in Discussion.
Table 1. WAME weights for the Atrium dataset. Weights in percent from estimate of differential expression using WAME on the Atrium dataset.
Table 2. WAME weights for the COPD dataset. Weights in percent from estimate of differential expression using WAME on the COPD dataset.
The analysis methods described in Background are applied to the data and pvalues and ranks computed. The respective probability plots are shown in Figure 3, demonstrating that there are substantial differences in the distribution of pvalues between the different statistics. Since correlations and unequal variances are observed, the model assumptions of the alternative standard methods do not seem to hold. The pvalues could thereby have become optimistic. On the other hand, it cannot be ruled out that the temporary assumption in WAME of no regulated genes makes its pvalues conservative, which could also partly explain the differences. These problems are studied below by use of resampled data.
Figure 3. Observed probability plots. Empirical distribution of pvalues compared to the distribution expected for nondifferentially expressed genes. The OLM and LIMMA curves largely coincide, as does the identity line and the WAME curve.
A common alternative to using the pvalues as measures of significance is to consider the ranking of the genes, induced by the pvalues or test statistics, and to select a fixed number of top ranked genes for further investigations. In Table 3 and Table 4 the concordance of the ranked lists are shown. The results from the included methods differ, for instance those from WAME compared to the other methods. This is not surprising since high correlations and highly unequal variances were identified by WAME, giving rise to highly unequal weights.
Table 3. Concordance of top lists in the Atrium dataset. Number of mutually included genes in the top100 lists in the Atrium dataset as determined by the different methods.
Table 4. Concordance of top lists in the COPD dataset. Number of mutually included genes in the top100 lists in the COPD dataset as determined by the different methods.
Resampled data
To examine closer the effect of violated assumptions of independence and identical distribution, we repeatedly selected two random subgroups of four arrays from within one group in the original data and performed tests between those groups. This was performed 100 times for the largest group in each of the two real datasets. Differentially expressed genes have unequal expected values in the two populations being sampled (cf. (2)). Since we now sample twice from the same condition, no differentially expressed genes exist.
Figure 4 shows the empirical pvalue distributions for the resampled COPD data analysed with the four methods, together with the respective average empirical distribution,
Figure 4. Probability plots. Empirical distributions of pvalues for LIMMA, weighted LIMMA, OLM and WAME from tests on 100 resamples from the COPD dataset. Average empirical distribution indicated. Since no signal is added, the curves should ideally follow the diagonal.
where F_{i }denotes the empirical CDF from the ith of the 100 resamples. For WAME, the pvalue distributions are very close to the expected uniform. For OLM, LIMMA and weighted LIMMA there is a high variability between the pvalue distributions and they are in many cases substantially different from the expected uniform. For WAME, OLM and LIMMA, the respective average empirical distribution is approximately correct, while for weighted LIMMA it is clearly optimistic. The results for the Atrium dataset (see Additional file 3) are very similar.
Additional file 3. Probability plots for the Atrium dataset. Empirical distributions of pvalues for LIMMA, weighted LIMMA, OLM and WAME from tests on 100 resamples from the Atrium dataset. Average empirical distribution indicated. Since no signal is added, the curves should ideally follow the diagonal.
Format: PDF Size: 487KB Download file
This file can be viewed with: Adobe Acrobat Reader
Evaluation of power
To evaluate the power of the tests in the studied datasets, a known regulation is added to randomly selected genes in one of the resampled groups, created according to the previous section. Thus, the noise is obtained from the real data and only the signal is synthetic. Ideally, the power can then be estimated by the proportion of differentially expressed genes that have a computed pvalue less than a fixed level. However, valid pvalues of the test statistics cannot be obtained from the respective models since, as demonstrated above, the corresponding assumptions are typically not valid. Ideally, the pvalues would be determined by the true null distribution of the respective test statistics, given the arraywise quality deviations. In the simulation study, the critical value of the test statistics are therefore estimated from the empirical distribution of the test statistic for the unregulated genes. This is used to estimate the power of the different statistics (details are given in Methods).
The power estimates for the different methods are shown in Figure 5, for a level 0.1% test. The 0.1% level yields approximately 22 false positives if relatively few genes are in fact differentially expressed. For WAME, Σ is estimated both before and after adding a signal to 2228 genes (10%), thereby substantially affecting the estimate of Σ (cf. Figure 6). The powers of the two versions are nevertheless very similar (difference less than 0.003) and only the latter version is included in the plot.
Figure 5. Estimated power. Estimated power in the simulated data for level 0.1% tests, based on resamples from the respective larger group in the Atrium and COPD datasets. Power is estimated at the marked points and spline interpolation is used in between.
Figure 6. Average empirical pvalue distribution for WAME under regulation. Average empirical pvalue distribution of the unregulated genes, calculated using WAME, when 0%, 0.1%, 1%, 5% and 10% of the genes have a log_{2 }differential expression of 1, i.e. a twofold change. When genes are regulated the estimate of Σ is biased, leading to conservative, nondiagonal curves.
When the covariancestructure matrix Σ is estimated in WAME it is assumed that no genes are differentially expressed. Figure 6 includes the average empirical distribution for the pvalues of the unregulated genes when different proportions of the genes have a log_{2 }differential expression of 1. It is clear that the distributions are biased for high proportions, giving conservative pvalues, which should be an effect of biased estimates of Σ.
The results from the studied datasets indicate (i) that WAME offers a relevant power increase compared to the included alternatives, (ii) that weighted LIMMA does not offer an advantage compared to LIMMA and (iii) that the moderated statistics (WAME, LIMMA and wLIMMA) are superior to the traditional methods of ranking by ordinary tstatistic (OLM) or estimated differential expression (FC).
Discussion
The WAME model and the simulations
The WAME model aims at catching quality deviations by one covariancestructure matrix common for all genes. This is certainly simplistic in some cases, e.g. when only certain physical parts of an array or certain types of mRNAs are of decreased quality. The estimated covariance structure can then only be expected to reflect a mixture of the qualities of the different genes. However, examining the simulations (Figure 5), we see a clear power gain in the WAME model compared to the other models. Also, WAME succeeds in catching enough of the quality deviations to make the pvalue distributions more correct, thus providing increased usefulness of the pvalues (Figure 3).
The models of LIMMA, weighted LIMMA and WAME are nested, where weighted LIMMA adds unequal variances and WAME adds unequal variances and correlations. Examination of Figure 1 shows that there are evident differences in variability between arrays. It is therefore interesting that we have not found a power increase of weighted LIMMA compared to LIMMA. Further, the pvalues of weighted LIMMA turned out to be too optimistic (Figure 4). Comparison with the results of the WAME method, where the power increases and the pvalue distributions get substantially more correct, suggests that the correlations are crucial in the model.
In the simulations, noise is taken from real data through resampling within a fixed group. This procedure provides data with fewer assumption on the noise structure compared to a fully parameterised simulation and should hopefully better reflect realistic situations. To evaluate the power of the different methods, a synthetic signal which is constant within each condition is added to the resamplebased noise. This follows the assumption in the models of both WAME, OLM, LIMMA and weighted LIMMA, that the noise structure is equal for genes that are differentially expressed and nondifferentially expressed. However, the biological variability of the expression of differentially expressed genes might be different under the different conditions due to the changed rôle of those genes. For complicated conditions such as complex diseases, the problem is more severe (cf. [1517]) since crucial genes might only be differentially expressed in a subset of the studied arrays. Further work is needed to evaluate the performance of WAME in such settings, as well as to possibly expand it to better fit these situations.
A relevant question regarding the modelling of quality deviations by the covariancestructure matrix Σ is whether biologically interesting features may be hidden by this model. In the present datasets, the question can partly be answered by examining the pairwise plots (cf. Figure 2) and noticing that a large proportion of the genes show similar deviations, which should speak against a specific interesting biological explanation. Also, the estimated covariance structure matrix Σ can be inspected with the goal of finding relevant correlations between arrays and thus highlighting interesting features in the data. Possible future work is to use such inspections to reveal unwanted features in normalisation or in preprocessing wetlab steps that give rise to correlated errors for a large proportion of the genes.
Weights with switched signs
In the studied datasets, strong correlations combined with unequal variances make some weights within a group switch sign, in essence meaning that it is beneficial to partly subtract some arrays within a group in the estimate to be able to add more of the others in the same group (cf. Table 1 and Table 2). Since this might seem counterintuitive, an elucidating example of possible mechanisms behind such weights follows.
Consider an example where two twocolour arrays are observed, X_{1 }and X_{2}. Let the two arrays have two sources of variation, one that is mutually independent (ε_{1}, ε_{2}) and one consisting of different proportions, a_{1 }and a_{2}, of one common source of variation η. Let ε_{1}, ε_{2 }and η be independent and normally distributed with expectation 0 and variances and , respectively. Furthermore, let μ be the parameter to be estimated. The model becomes
Then, X_{1 }gets a negative weight if and only if
i.e. if array 1 includes a large enough contribution from the common source of variation. When a negative weight is allowed instead of removing the array, a smaller proportion of the common source of variation is included in the final estimate. Its precision is thus increased.
Validity of the pvalues and derived entities
Varying quality of arrays and correlated errors were demonstrated in [2,3] and in the present paper through examination of the data. These questions are typically neglected in microarray analyses, both when using parametric and when using nonparametric analysis methods, since independence and identical distribution or exchangeability are generally assumed under the null hypothesis. Thus, the validity is questionable of the corresponding pvalues and their derived entities, e.g. false discovery rates and estimates of proportions of differentially expressed genes. This problem is obvious in the resample based simulations.
A number of experiments have been analysed (data not shown) in addition to those published in the present paper and in [2,3]. In almost all cases relevant unequal variances and correlations have been identified, indicating that the problem is common.
In the resample based simulations with added signal, WAME is shown to be conservative, which is an effect of the biased estimate of Σ. Further work on an estimator of Σ with better characteristics under regulation is therefore needed. However, the simulations indicate (i) that the power of the test is basically unaffected by the bias and (ii) that hundreds of genes may be differentially expressed (twofold) with only mildly conservative pvalues as result.
Correlations between genes or between arrays?
It has recently been argued that the expression of different genes are highly dependent, making the law of large number normally inapplicable [18] and standard estimators of e.g. the false discovery rate (FDR) imprecise [19]. In [19], a latent FDR is introduced, being the conditional FDR given a random effect b that captures the correlation effects between genes. The FDR is then the marginal latent FDR, that is the average over the random effect b.
For the datasets examined in the present paper, the model assumptions of e.g. the ordinary linear model are shown not to hold (cf. Figure 1 and Figure 2). This can be expected to result in invalid pvalues, which is indeed observed in Figure 4. Interestingly, the pvalue distribution seem to be valid marginally, i.e. on average over the different resamples, which would yield valid but imprecise estimates of the FDR. This type of failed model assumptions is not taken into account in e.g. [18,19]. Since for a performed experiment, the pvalues from the ordinary tstatistic (OLM) share a common bias conditional on the experiment (see Figure 4), the different pvalues may be highly dependent. However, this dependency is due to failure of taking arraywide quality deviations into account in the model and not due to the nature of microarray data per se, e.g. through substantial longrange genegene interactions.
Consequently, the strong observed dependencies between statistics from different genes might largely be explainable by quality deviations between the arrays in the experiment, e.g. correlations between arrays. Since WAME models these deviations such that the pvalues are essentially correctly distributed when few genes are differentially expressed in the studied datasets, the dependency between genes should be greatly decreased. The covariance structure matrix Σ is therefore in a sense a parallel to the random factor b in [19]. It remains as future work to evaluate the genegene dependencies and estimates of e.g. the FDR in the context of the WAME model.
In the WAME model, the data from different genes are assumed independent, which is unrealistic, e.g. since genes act together in pathways. However, this is only used in the derivation of the maximum likelihood estimaties of the covariance structure matrix Σ and the shape parameter α. The assumption could thus be relaxed to a dependence between the different genes that is weak enough that the estimates of Σ and α become precise, and accurate under H_{0}. This holds if the law of large numbers is applicable for averages of certain functions of the genewise observed data (cf. the likelihood functions in [2,3]). Given the large number of genes and the observed pvalue distributions in Figure 4, this relaxed assumption seems plausible.
It can be noted that for the studied data, WAME has higher power and considerably more valid pvalues than weighted LIMMA. Since the difference between the weighted LIMMA and WAME models is the inclusion of correlations between arrays, this emphasises the importance of the correlations in the model.
Conclusion
Statistical methods in microarray analysis are typically based on the often erroneous assumption that the data from different arrays are independent and identically distributed. An exception is Weighted Analysis of Microarray Experiment (WAME) where heteroscedasticity and correlations between arrays are modelled by a covariancestructure common for all genes. In the present paper, WAME has been extended to handle datasets without a natural pairing, e.g. data from onechannel microarrays, and corresponding estimates and test statistics have been derived. In the examined onechannel microarray datasets WAME detected unequal variances and nonzero correlations.
WAME was compared with four other common methods: an ordinary linear model with ttests, LIMMA, weighted LIMMA, and foldchange ranking. The comparison was performed using resampling of the different arrays within the datasets. Here, WAME had the highest power. When a relatively small proportion of the genes are regulated, WAME also generates close to correct pvalue distributions while the pvalue distributions from the other methods are highly variable. However, when many genes are differentially expressed, the pvalues from WAME tend to be conservative.
In conclusion, pvalues from the standard methods for microarray analysis should in general not be trusted and any result based on pvalues, e.g. estimates of the number of regulated genes and false discovery rates, should be interpreted with care. The analyses of the examined datasets showed that WAME gives a powerful approach for finding differentially expressed genes and that it produces more trustworthy pvalues when a relatively small proportion of genes are differentially expressed.
Methods
Details on the new version of WAME
Model Framework
For g = 1,..., m, let X_{g }be an ndimensional vector with expectation μ_{g }= Dγ_{g}, where D is the design matrix, having rank k, and γ_{g }∈ ℝ^{q }is the parameter vector. Furthermore, let
where Σ is the nonsingular covariancestructure matrix, c_{g }is the variancescaling factor, α is the shape parameter for c_{g }and (c_{1}, X_{1}),...,(c_{m}, X_{m}) are assumed independent. The differential expression vector is defined as
where C is a matrix of rank p such that δ_{g }is estimable. Here, an estimator of δ_{g }and a test for
are in focus.
As mentioned in Background, one main obstacle is that Σ is hard to estimate. In fact, Σ and δ_{g }cannot be maximum likelihood estimated simultaneously, since there are trivial infinite suprema of the likelihood, e.g. when the variance of one observation is set to zero and the corresponding mean is selected so that it equals that observation.
The current WAME method
In the current version of WAME [3], Σ is estimated as follows. First, temporarily assume that μ_{g }= 0 for all genes, which is reasonable for paired experimental designs with few differentially expressed genes between any pairwise measured conditions. For each gene, the variance scaling factor c_{g }is removed by dividing the n measurements with the first measurement, yielding a random vector distributed according to a multivariate generalisation of the Cauchy distribution. A scaled version of Σ is then maximum likelihood estimated numerically. Second, the unknown scale and the hyperparameter α of the prior distribution of c_{g }are maximum likelihood estimated numerically without the assumption of μ_{g }= 0. The parameters Σ and α are subsequently treated as known in the maximumlikelihood estimates and likelihoodratio tests for the different genes.
The new WAME method
The new version of WAME relaxes the assumption from μ_{g }= 0 to δ_{g }= 0, which incorporates many designs without a natural pairing. This is performed by subtracting an arbitrary estimator of μ_{g}, which is unbiased under H_{0 }and has as image the space _{0 }of possible values for μ_{g }under H_{0},
It can be shown that this transformation preserves the estimability of δ_{g}.
By construction, the transformed data Y_{g }will have expectation zero for nondifferentially expressed genes and the current WAME method can be applied on Y_{g}, including the estimation of the covariancestructure matrix Σ_{Y }for Y_{g}. It will now be proved that the likelihood ratio tests of (9) and the maximum likelihood estimates of δ_{g }based on X_{g }or Y_{g }are identical, if α and Σ or Σ_{Y }respectively are considered known.
We shall henceforth consider a fixed gene g and drop the g index.
Equality of tests and estimators
Before beginning, some further definitions are needed. Define the Mahalanobis inner product corresponding to a symmetric n by n matrix A as
and the norm ∥·∥_{A }as
where x, x_{1}, x_{2 }lies in the rowspace of A and the generalised inverse A^{}is any matrix satisfying AA^{}A = A. Let χ denote the ndimensional inner product space with ⟨·,·⟩_{Σ }as inner product. Define ⊂ χ as the space of possible values for μ_{g},
and let _{0 }⊂ χ denote the corresponding space restricted by the null hypothesis,
Proposition Let be an arbitrary linear estimator of μ, which is unbiased under H_{0 }and which has image _{0}. Let
and let Σ_{Y }be the covariancestructure matrix of Y. Then the likelihood ratio test of (9) and the maximum likelihood estimate of δ based on X with Σ and α known are identical to the ones based on Y with Σ_{Y }and α known.
Proof of the Proposition
The proof is divided into two steps which combined conclude the proof.
1. The likelihood ratio test (LRT) of (9) and the maximum likelihood estimator (MLE) of δ does not depend on the choice of .
2. The proposition holds when is the MLE of μ under H_{0}.
Proof of step 1
Let μ' and μ'' be two valid choices of , i.e. they are both unbiased estimators of μ under H_{0 }and have _{0 }as image. Let Y' = X  μ' and Y'' = X  μ''. Recall that a matrix P is a projection matrix projecting on _{0 }if and only if for all x ∈ ℝ^{n}, Px ∈ _{0 }and for all x_{0 }∈ _{0}, Px_{0 }= x_{0}. It can be shown that μ' and μ'' can be written as μ' = P' X and μ'' = P'' X for some projection matrices P' and P'' projecting on _{0}. Since P' and P'' project on the same space it follows that P' P'' = P'' and P'' P' = P', and thus (I  P') Y'' = Y' and (I  P'') Y' = Y''. Hence there is an invertible map between Y' and Y'' and likelihood methods based on Y' and Y'' respectively will give equal results. Consequently, the MLE of (9) and the LRT of δ will not depend on the choice of
Proof of step 2
Since δ is estimable based on X, there exist a matrix A such that C = AD and thus δ = Aμ. The likelihood of μ can therefore be examined instead of the likelihood of δ.
The likelihood of μ based on X can be shown to be
where ∝ denotes proportionality. Using the Projection Theorem [20], the MLE of μ is the orthogonal projection of X on ,
where the orthogonality is according to the inner product of χ. When H_{0 }is true, μ is restricted to _{0 }and thus the MLE of μ becomes
Note that is a valid choice for , i.e. is unbiased under H_{0 }and has _{0 }as image. Let
which gives , where denotes the orthogonal complement of _{0 }in χ. Standard properties of the normal distribution gives
where μ_{z }= D_{z}γ with D_{z }= , and where .
The likelihood function of μ_{z }(with respect to the Lebesgue measure on the space of possible values of Z spanned by the column vectors of Σ_{z}) can be written as
Since, δ is estimable based on Z, the likelihood of μ_{z }can be examined instead of the likelihood of δ.
The likelihood ratio statistic of (9) based on X is defined by
which can be rewritten (cf. [3]) as a strictly increasing function of
where ^{⊥ }and are the orthogonal complements of and _{0 }respectively.
Note that the space of possible values for μ_{z }is ∩ and that μ_{z }= 0 under H_{0}. Let denote the orthogonal projection according to . Then, the likelihood ratio statistic of (9) based on Z can in analogy with (13) be shown to be a strictly increasing function of
The Lemma below yields that for all and all z ∈ , and . The equivalence of the test statistics (13) and (14) is now straightforward,
Lemma Let be a subspace of χ and let be the orthogonal projection from χ onto . Then for any x_{1}, x_{2 }∈ ,
Proof Let A be a matrix of a change of basis [20] from the standard basis to an orthonormal basis of χ such that the first l basis vectors span . Let I_{(l) }denote the identity matrix with all but the l top left diagonal elements set to zero. It follows that A^{T }A = Σ^{1 }and A = I_{(l)}A and therefore,
where the last equality uses the fact that (AB)^{ }= B^{ }A^{1 }when A is invertible □.
The next step is to show that the MLE of δ when X is observed is identical to the MLE of δ when Z is observed. The former is defined by
Define _{0 }= {γ : Dγ ∈ _{0}} and _{1 }= {γ : Dγ ∈ } and note that for any γ there exist γ_{0 }∈ _{0 }and γ_{1 }∈ _{1 }such that γ = γ_{0 }+ γ_{1}. Thus,
where the second equality follows from the generalised Theorem of Pythagoras [20], the Lemma, and the fact that and . Now since γ_{0 }and γ_{1 }minimise the expression independently of each other and since Cγ_{0 }= 0 by construction,
For all γ_{0 }∈ _{0}, Cγ_{0 }= 0 and D_{z}γ_{0 }= 0, so the area of minimisation can be extended,
which is the MLE of δ based on Z by definition. □
Remark 1 Using the invertible map between any two choices of Y, Y and Y', as defined in Step 1 above, the respective maximum likelihood estimates of α, Σ_{y }and Σ_{y' }can be shown to be identical based on Y or Y'. In this sense, the choice of is thus truly irrelevant.
Remark 2 Sometimes, additional linear combinations of γ can be assumed to be zero for most genes, C* γ = 0 for some matrix C* with rowspace being a superspace of the rowspace of C. Let P* be any projection matrix on the corresponding space _{* }= {μ : μ = Dγ, C* γ = 0, γ ∈ ℝ^{q}} and let Y* = X  P* X. It is straight forward to show that a variant of the Proposition still holds, so given the covariance structure matrices the inference results concerning Cγ will be identical, based on Y or Y* respectively. However, the estimates of the covariance structure matrices for Y and Y* might not be coherent and the results are expected to differ slightly.
The estimator of power
Consider a certain experimental design, a level 1λ test and a differential expression δ. Let a realisation of the experiment be given, which e.g. results in certain quality deviations between arrays. The conditional power is defined as the probability of identifying a random gene in the current experiment, i.e. conditional on e.g. the quality deviations, when the gene has differential expression δ. The power is then defined as the average conditional power over all possible realisations of the experimental design. The power is thus related to an unperformed experiment while the conditional power is related to a specific performed experiment. Here, the test is assumed to be valid conditional on the experiment, i.e. to have conditional power λ when δ = 0.
In Evaluation of power, the aim is to estimate the power for a hypothetical experiment where the distribution of the data under the null hypothesis is obtained by resampling of real data. For a given resample, a constant differential expression is added to randomly selected genes and the statistics t_{g }are computed. The estimate of the conditional critical value is computed so that a proportion λ of the unregulated genes satisfy t_{g} ≥ . The conditional power is then estimated by the proportion of regulated genes satisfying t_{g} ≥ . The power is finally estimated by averaging the estimated conditional power over the resamples.
Competing interests
The author(s) declares that there are no competing interests.
Authors' contributions
ON provided initial ideas related to the generalisation. AS formulated the generalisation, performed the proofs in Methods, designed and programmed the analyses, simulations and plots. EK and ON helped refining the generalisation and Methods. AS, EK and MR drafted the manuscript. All authors continuously provided feedback on various parts of the work leading to the manuscript and approved the final version of the manuscript.
Acknowledgements
Lina Gunnarsson is acknowledged for providing feedback on the manuscript. AS and EK acknowledge the National Research School in Genomics and Bioinformatics for funding.
References

Auer H, Lyianarachchi S, Newsom D, Klisovic M, Marcucci G, Kornacker K: Chipping away at the chip bias: RNA degradation in microarray analysis.
Nature Genetics 2003, 35:292293. PubMed Abstract  Publisher Full Text

Kristiansson E, Sjögren A, Rudemo M, Nerman O: Weighted Analysis of Paired Microarray Experiments.
Stat Appl Genet Mol Biol 2005, 4:Article30. PubMed Abstract  Publisher Full Text

Kristiansson E, Sjögren A, Rudemo M, Nerman O: Quality Optimised Analysis of General Paired Microarray Experiments.
Stat Appl Genet Mol Biol 2006, 5:Article10. PubMed Abstract  Publisher Full Text

Arnold S: The Theory of Linear Models and Multivariate Analysis. John Wiley & Sons; 1980.

Smyth G: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
Stat Appl Genet Mol Biol 2004, 3:Article3. PubMed Abstract

Ritchie M, Diyagama D, Neilson J, van Laar R, Dobrovic A, Holloway A, Smyth G: Empirical array quality weights in the analysis of microarray data.
BMC Bioinformatics 2006, 7:261. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

R Development Core Team: [http://www.Rproject.org] webcite
R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria; 2006.
[ISBN 3900051070]

The WAME homepage [http://wame.math.chalmers.se] webcite

Edgar R, Domrachev M, Lash A: Gene Expression Omnibus: NCBI gene expression and hybridization array data repository.
Nucleic Acids Research 2002, 30:207210. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Barth A, Merk S, Arnoldi E, Zwermann L, Kloos P, Gebauer M, Steinmeyer K, Bleich M, Kaab S, Hinterseer M, Kartmann H, Kreuzer E, Dugas M, Steinbeck G, Nabauer M: Reprogramming of the Human Atrial Transcriptome in Permanent Atrial Fibrillation.
Circulation Research 2005, 96(9):10221029. PubMed Abstract  Publisher Full Text

Spira A, Beane J, PintoPlata V, Kadar A, Liu G, Shah V, Celli B, Brody J: Gene expression profiling of human lung tissue from smokers with severe emphysema.
Am J Respir Cell Mol Biol 2004, 31(6):60110. PubMed Abstract  Publisher Full Text

Irizarry R, Hobbs B, Collin F, BeazerBarclay Y, Antonellis K, Scherf U, Speed T: Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data.
Biostatistics 2003, 4(2):24964. PubMed Abstract  Publisher Full Text

Gentleman R, Carey V, Bates D, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Friedrich L, Li C, Maechler M, Rossini A, Sawitzki G, Smith C, Smyth G, Tierney L, Yang J, Zhang J: Bioconductor: Open software development for computational biology and bioinformatics.
Genome Biology 2004, 5:R80. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Tomlins S, Rhodes D, Perner S, Dhanasekaran S, Mehra R, Sun X, Varambally S, Cao X, Tchinda J, Kuefer R, Lee C, Montie J, Shah R, Pienta K, Rubin M, Chinnaiyan A: Recurrent Fusion of TMPRSS2 and ETS Transcription Factor Genes in Prostate Cancer.
Science 2005, 310:644648. PubMed Abstract  Publisher Full Text

van Wieringen W, van de Wiel M, van der Vaart A: A Test for Partial Differential Expression. In Tech Rep WS20064. Department of Mathematics, Vrije Universiteit; 2006.

Tibshirani R, Hastie T: Outlier sums for differential gene expression analysis.
Biostatistics 2007, 8(1):28. PubMed Abstract  Publisher Full Text

Klebanov L, Yakovlev A: Treating Expression Levels of Different Genes as a Sample in Microarray Data Analysis: Is it Worth a Risk?
Stat Appl Genet Mol Biol 2006, 5:Article9. PubMed Abstract

Pawitan Y, Calza S, Ploner A: Estimation of the false discovery proportion under general dependence.
Bioinformatics 2006, 22(24):30253031. PubMed Abstract  Publisher Full Text

Anton H: Elementary Linear Algebra. 6th edition. Wiley; 1991.