Statistical methodology for the analysis of dye-switch microarray experiments

Background In individually dye-balanced microarray designs, each biological sample is hybridized on two different slides, once with Cy3 and once with Cy5. While this strategy ensures an automatic correction of the gene-specific labelling bias, it also induces dependencies between log-ratio measurements that must be taken into account in the statistical analysis. Results We present two original statistical procedures for the statistical analysis of individually balanced designs. These procedures are compared with the usual ML and REML mixed model procedures proposed in most statistical toolboxes, on both simulated and real data. Conclusion The UP procedure we propose as an alternative to usual mixed model procedures is more efficient and significantly faster to compute. This result provides some useful guidelines for the analysis of complex designs.


Background
DNA microarray technology is a high throughput technique by which the expression of the whole genome is studied in a single experiment. Experiments must be well organized and design issues are crucial, see [1,2]. In dual label experiments Cy3 and Cy5 are used as fluorescent dyes allowing to compare two RNA samples on the same slide. It is now well known that there exists a differential effect of the two dyes [3,4], that can be gene-specific. An efficient way to remove this technical artifact is to use balanced reverse dye designs [5]. Balanced reverse dye designs can be divided into three classes along a line of strengthening balancing constraints: 1. Balanced reverse dyes for which each biological sample is hybridized only one time and therefore present with only one dye, on only one array (Table 1.1). These designs are globally balanced but not individually balanced.
2. Individually-balanced design for which each biological sample is divided into two parts, one hybridized with Cy3 on one array and the other with Cy5 on another array. Each biological sample is hybridized exactly two times (Table 1.2). 3. Dye-swaps for which each couple of biological samples from two conditions are hybridized on two arrays with reverse dyes. Dye-swaps are constrained to be couple-balanced (Table 1.3).
Dye-swap design is mostly used when the technical error is higher than the biological variability, either to reduce the technical variance, or when gene-specific dye-bias is of concern [6,7]. When the biological variability is greater than the technical error, globally balanced designs are statistically more efficient [5]. However the number of biological samples is sometimes limited, therefore this design is not always possible in practice.
The term Dye-switch is used for the first and sometimes also for the second classes. Dye-switch designs of the second class are sometimes described and proposed in papers dealing with microarrays experiments. For example loop designs are often members of this class [8,9], although the distinction between the first and the second class is not always clearly made.
A major point to notice is that the statistical analysis may be very different for the three classes of design. The analysis of the first and third classes is straightforward and well described in articles (see for example [4,10,11]): the experimental units are mutually independent (we consider as usual that the two conjugate arrays of the dyeswaps are summed up to only one experimental unit), and simple statistical procedures such as Student T-tests (or regularized T-tests) can be performed. On the contrary, if we consider the second class of designs, the experimental units are not independent, a feature that must (or must not) be accounted for. The literature about the statistical study of such designs is limited: some papers proposed some theoretical contributions for their analysis [12,13], but simple guidelines for experimenters and practical considerations (computational burden, choice of a strategy for parameter estimation) are not available.
We consider here the simplest individually-balanced dyeswitch design: two conditions A and B are compared in a two-color cDNA microarray experiment, with n biological samples for each condition. The design is the following: each RNA sample (A 1 to A n for condition A, and B 1 to B n for condition B) is divided into two parts, one labelled with Cy5 and the second labelled with Cy3. Then 2n microarrays are hybridized with respectively A 1 Cy5 and B 1 Cy3, B 1 Cy5 and A 2 Cy3, A 2 Cy5 and B 2 Cy3, and so on till B n Cy5 and A 1 Cy3, (see Table 1.2). There are 2n samples, 4n labelled samples, 2n microarrays, and each sample is hybridized two times (one with Cy5 and one with Cy3) on two different arrays. We propose a simple, efficient and robust method for the statistical analysis of this experiment.

Model on the measure of the expression of genes
After the normalization step, X i is the expression measure on the log-scale, for a given gene, corresponding to condition A on array i. Let j(i) denote the sample number corresponding to condition A and array i.
Similarly, Y i is the expression measure for the condition B sample on the same array, and j'(i) the sample number corresponding to condition B and array i. In the following the gene index is not present in order to simplify the mathematical expressions, but it is important to note that all the terms in the following models are gene-specific.
Here we use an analysis of variance (ANOVA) model for the expression measure as introduced by [10].
The model for X i and Y i is the following: where • µ A and µ B are the population mean expression measures for condition A and B.
• δ l(i) is a two-level fixed effect corresponding to the dye effect. δ l(i) = δ 1 (resp. δ 2 ) for all the samples labelled with Cy5 (resp. Cy3). This term accounts for the gene-specific dye bias.
• B j(i) represents an independent gaussian random term with mean 0 and standard deviation σ B , corresponding to the random effect of sample j(i). This variable is specific to the biological sample and is called biological error, related to the variability of the biological material inside each population A and B.
• M i represents an independent gaussian random term with mean 0 and standard deviation σ M . M i is the effect of the spot associated to the gene under concern in microarray i and has the same value for the two samples which are hybridized on array i. This error term takes into account the spatial heterogeneity in each array that affects both Cy3 and Cy5 measurements.
• T i represents an independent gaussian random term with mean 0 and standard deviation σ T , corresponding to the technical variability, including the steps of labelling, hybridization and measure of intensity of fluorescence. This variable has a specific value for each combination gene×dye×sample, even if the samples are hybridized on the same array and at the same spot, so that T i and are independent random variables. T i and M i are the two components of the so-called technical error.

Model on the difference of expression on one array
..,2n. Using equation (1) we obtain: which may be written is the true differential expression between conditions A and B for the gene under concern, • BD i = B j(i) -B j'(i) is a random variable with mean 0 and standard deviation σ B , • TD i = T i -is an independent random variable with mean 0 and standard deviation σ T , • δ = δ 1 -δ 2 is the difference between the Cy3 and Cy5 dye effects. This term accounts for the gene-specific dye bias.
Each variable D i follows a Gaussian distribution with mean E(D i ) = µ + (-1) i+1 δ and variance . All the covariances cov(D i , D j ) are equal to zero except the following ones: with the convention that 2n + 1 = 1.
In this study, we present and compare different strategies for the statistical analysis of individually-balanced designs. The article is organized as follows. In the Results section, five statistical procedures to analyze individually balanced designs (Table 1.2) are compared on both simulated and real data. The Conclusion section draws the main conclusions and gives some useful guidelines for the analysis of individually-balanced designs. The details of the computation are given in the Methods section.

Statistical procedure comparison
In this section, we investigate the efficiency of five test procedures for the differential analysis of datasets corresponding to the design of Table 1.2. The procedures are the following (see the Methods section for more details): • Naive Method NM: for each gene, the naive test statistic is computed.
• Unbiased Paired Method (UP): for each gene, the unbiased paired statistic is computed. Notice that from the Methods section, the theoretical value of C must be positive. In practice, the estimated value may be negative. In such a case, C is truncated at 0.
• Unbiased Unpaired Method (UU): for each gene, the unbiased unpaired statistic is computed. As for the previous method, the value of C XY must be positive. If not, C XY is truncated at 0. Furthermore, the unbiased variance estimator is . Since C XY is non-negative, the variance estimator may have a negative value. In such a case, the variance can be fixed at a given threshold (0.001 in the following). It is important to consider both the ML and REML algorithms for the mixed model since each algorithm has its own advantages. While ML is known to provide biased estimates of the variance components, computations are faster and the algorithm does converge. REML gives unbiased estimates of the parameters, but may not converge if the number of observations is small. Both ML and REML computations were performed using the R package Maanova [10].

Simulations
To study the behavior of the 5 procedures, we performed a simulation study using model (1)

Control of the Type I error rate
We first consider the case µ = 0. Table 2 shows the actual Type I error rate level of the 5 test procedures, when the requested nominal level is 5%. Different behaviors can be observed: NM and ML result in a type I error rate higher than the nominal level, and procedure UU is conservative. UP results in an actual level that is close to the expected one, whatever the conditions. In most cases, REML enables an efficient control of the type I error. Yet, when the biological variability is high and the number of samples is low, REML yields a high type I error because of inconsistent estimations of the variance (see the next section). When = 2 and n = 5, the discrepancy between the theoretical and the actual level is even worse for REML than for the other methods.
From these first observations we conclude that we can discard procedures NM and ML, since in differential analysis an effective control of the Type I error rate is necessary.

Performance analysis
We now compare the performance of the 3 remaining procedures to detect differentially expressed genes. Table 3 shows the proportion of detected differentially expressed genes, for different values of the parameter set. It clearly appears that the power of procedure UU is low compared with procedures UP and REML. This may be the consequence of the Student approximation (each test statistic is compared with the quantile of a Student distribution with 2n -2 degrees of freedom), that could be more erroneous in the case of the UU statistic.
An interesting point is that UP results are more stable than the REML results. If we consider sample sizes n larger than 20, we observe that the absolute values of the approximate REML T-test range from 0 to 32, except for some genes where the absolute value is larger than 400. These outliers come from an erroneous estimation of the variance of the mean difference, that is evaluated to be (almost) 0. This does not happen with (UP) since the estimated variance is max(S 2 , S 2 + 2C), i.e. the variance is overestimated to avoid outliers. Notice that despite this overestimation in many cases the power of UP is larger than the power of REML.

Computational burden and convergence
We now consider the important question of computational time for the 2 competitive procedures UP and REML. Since microarray experiments can involve hun-   Table 4 gives the user CPU time associated to each procedure for the complete analysis of 10,000 genes. While the computational time is constant whatever the condition for the (UP) procedure, (REML) is 8 to 330 times longer than (UP), depending on the number of samples.
Furthermore, REML can result in inconsistent estimates of the variance, as shown in the previous sections, or may not converge. Table 4 provides the number of genes for which the REML algorithm did not converge.

Embriogenomic data
The impact of pregnancy on the cattle endometrium transcriptom is investigated in [14]. In Mammals, the implantation of the embryo is a key event in the establishment of a pregnancy. A microarray experiment has been made to analyze the gene expression of the bovine pregnant endometrium and determine key pathways that control the endometrium physiology during the implantation process. The expression of 13300 genes in the endometrium of cows (n = 5) has been investigated. Only 5 animals were available for each condition so that the dye-switch design of Table 1.2 was used. Gene profiling has been established to analyze the impact of pregnancy by comparing the endometrium of cyclic (day 20 of cycle) versus pregnant animals (day 20 of pregnancy). In the following, the results of the five statistical procedures defined above are compared using this dataset.
The Venn diagram of Figure 1 shows the number of genes declared differentially expressed (DE) by 4 methods using the Bonferroni method with a 5% level. The UU method gives the least number of DE genes (4) and is not presented in the diagram. REML (which did not converge for 3 genes) gives the greater number of DE genes (93), among which 23 are also found by the other methods, and 70 are specifically found by REML (70 REML specific genes). 70 genes are found DE by ML (22 ML specific genes), and 58 by the naive method (9 Naive specific).
Finally 33 genes are declared DE by UP, and all of them are also found by one, two or all of its competitors. Therefore UP provides the less discordant results. The higher number of DE genes obtained with the naive and the ML methods was expected, since it is known from the theory and the simulation study that these methods yield more false positives than the nominal risk. Figure 2 (right) shows that the ML and UP estimates of the standard error are coherent but that the ML estimate are lower than the ones obtained by the UP method. This point is in keeping with the statistical theory which assesses that the UP estimate of the variance is unbiased while the ML estimate has a negative bias.
The high number of DE genes specifically found with REML is odd. Figures 2 and 3 show that this comes from very low estimates of variance for some genes, so that these genes are declared DE not because the mean difference of expression between the two conditions is high, Venn diagram for the embriogenomics experiment   but because this mean difference is divided by a very low standard error. So most of the 70 genes only found by the REML method are due to too low estimates of the gene variance obtained by the REML algorithm. This observation is in keeping with the results of the simulation study. Therefore the UP method or the naive method should be preferred in this particular experiment. The use of REML without a sharp biological analysis of the results gene by gene would be misleading.

Teleost fish dataset
An important application of the methodology proposed in the previous section is the analysis of loop design experiments. Loop and interwoven loop designs were initially proposed in [2] to compare p treatments, where p is 3 or higher. Figure 4 displays a particular interwoven loop design where 3 different 2-by-2 loop comparisons of treatments are combined in a single experiment. The 3 loop comparisons are Each of these comparisons corresponds to the design of Table 1.2 discussed in the previous section. Such experimental designs have been studied both theoretically [15] and practically [8,9]. Here, we briefly present the Teleost fish data of [8].
The Teleost fish experiment aims to compare 3 populations of fish (Northern Fundulus heteroclitus, Southern Fundulus heteroclitus and Fundulus grandis). Five individuals were examined in each population to determine the variation in gene expression between populations. Each individual is used to probe four cDNA microarrays, according to the design of Figure 4. The raw data consist of 120 measurements (15 individuals × 4 slides × 2 duplicates per slide) for 907 genes.
In [8], the signal is modelled as follows (after per slide duplicate averaging): where A, D, G and V stand for Array, Dye, Gene and Variety, respectively. Then the 4 measurements corresponding to a given individual are averaged, and an F statistic is computed per gene to check whether the variety effect is significant or not.
This strategy roughly amounts to the UU test procedure of section when the number of treatments is higher than 2.
The main difference is that in model (4), the model does not include the array random effect which takes into account the dependency between two measures on the same array. According to the results of section, this implies that the variance estimator is biased, leading to a loss of power.
As an alternative, we perform the statistical analysis using the UP procedure. Each pairwise comparison between 2 varieties is made, and a gene is declared differentially expressed if at least 2 of the 3 tests are significant. Each test is performed at the level 0.02, meaning that for a given gene, the nominal level is roughly 0.001 (3 × 0.02 2 for 2 of the 3 tests to be significant under H 0 at level 0.02). This is a good compromise between the 0.01 threshold adopted in the original articles with no correction for multiple testing, and the 0.5 × 10 -4 (0.05/907) threshold given by a 5% level per test combined with a Bonferronni multiple testing correction. While the drawback of our strategy is to replace one test by three, the advantage is that the variance estimate is unbiased.
Comparison of the standard errors obtained with ML, REML and UP for the REML-DE genes of the embriogenomics experiment  Table 5 gives the Oleksiak original list of differentially expressed genes found with the original method and the UP list of genes found with the UP procedure.
Among the 15 genes originally identified, 5 are also declared differentially expressed with the (UP) method. At a first glance, the (UP) procedure seems less powerful than the original method since only 9 genes are found here compared with the 15 genes of the original article. But due to the threshold adopted by the authors in [8], the expected number of false positives is 9 for the Oleksiak list, whereas for the (UP) list we expect only 1 false positive. Therefore most of the 10 extra genes found in [8] may be false positives. To examine the discriminant effect of the 9 genes of the (UP) list, we performed as in the original publication a clustering of the individuals, according to the significative genes. The corresponding tree is given in Figure 5. With only 9 genes (rather than 15), the classification obtained with (UP) is improved compared with the classification of the original method.

Discussion
Random terms taking into account the array and the sample effects must be included in the statistical model at the gene level for dye-switch experiments. We showed on simulations that the naive paired T-test, which does not take into account the biological sample effect, leads to more false positives than expected, especially when the biological sample effect is high. It may be safely used only when the biological variance is lower than the technical variance. The REML estimate for mixed model provides an approximatively correct false positive rate, at the price of high computational complexity, lack of convergence for low or medium sample sizes and sometimes spurious results. To the contrary, the UP method we propose is easy to implement and not computationally intensive. The method is protected against spurious results, leading to a more robust and powerful analysis than REML when the biological variability is high and the number of samples low, an usual situation in microarray experiments.
Mean difference versus standard error for the REML-differentially expressed genes of the embriogenomics experiment Figure 3 Mean difference versus standard error for the REML-differentially expressed genes of the embriogenomics experiment. Standard error of the difference obtained by REML (y-axis) versus mean difference between the two conditions (x-axis). Black points are not found DE by other methods than REML.
For small sample size experiments, it is advised to use regularized T-test, see [16][17][18][19]. Regularization strategies are based on statistical methods that take the individual variance of each gene as input and give a regularized variance for each gene as output. The UP procedure proposed in this paper gives an estimate for the variance of the differential expression for each gene, so it allows a further application of all these regularization methods.

Conclusion
In this paper the proposed estimate of the variance of the differential expression is assessed for the comparison between two conditions in a dye-switch design. The same methodology could be extended to more complex designs involving more than two conditions and duplicate hybridizations of the same biological sample on different arrays.

Paired test procedure
According to expression (2), an unbiased estimator of µ is . The variance of this estimator is To perform a statistical test on parameter µ, we need to estimate .

Naive variance estimate
The The Teleofish experiment design. Lists of genes whose expression was significantly different between populations. The first 5 genes are found differentially expressed by both methods.
The expectation of this alternative estimator is S 2 is a downward biased estimator of V( ). The higher compared with , the higher the bias: From this naive estimate of the variance we can derive a first T-test statistic to be used for the differential analysis:  Clustering tree for the Teleofish dataset Figure 5 Clustering tree for the Teleofish dataset. Clustering tree for the Teleofish dataset, obtained from the second list of differentially expressed genes of Table 5.
which is approximately distributed as a Student distribution with 2n -2 df under H 0 .

Unpaired test procedure
Let (respectively ) be the mean of the 2 results obtained with the same biological sample (in 2 different arrays and with the 2 dyes) for condition A (respectively condition B). From model (1) one obtains where j is the biological sample index (recall that sample j is different for the two conditions), i(j) and i'(j) are the arrays on which sample j has been hybridized. and may be correlated as a result of a possible common array effect. and are uncorrelated because the two different biological samples of the same condition cannot be present together on the same array. From result (5) we have: The usual unpaired estimate of is equal to , where and , whose common mean (under the homoscedastic model (1)) is equal to Therefore This method overestimates the true variance The overestimation is more dramatic as increases.
This estimate may be corrected: may be estimated using the empirical covariance between and . Let with the convention that . The mean of the first sum is It is easy to see that the second sum in C XY has the same mean. Therefore an unbiased estimate of is , and the approximate unpaired t-statistic is

Authors' contributions
TMH and JJD conceived the method and prepared the manuscript. TMH and JA implemented part of the software and performed the statistical analysis. NM made the Embriogenomics experiment under the direction of OS. All authors contributed to the discussion and approved the final manuscript.