Statistical methodology for the analysis of dyeswitch microarray experiments
 Tristan MaryHuard^{1}Email author,
 Julie Aubert^{1},
 Nadera MansouriAttia^{2},
 Olivier Sandra^{2} and
 JeanJacques Daudin^{1}
https://doi.org/10.1186/14712105998
© MaryHuard et al; licensee BioMed Central Ltd. 2008
Received: 25 June 2007
Accepted: 13 February 2008
Published: 13 February 2008
Abstract
Background
In individually dyebalanced 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 genespecific labelling bias, it also induces dependencies between logratio 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.
Keywords
Background
 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.
Individuallybalanced design for which each biological sample is divided into two parts, one hybridized with Cy 3 on one array and the other with Cy 5 on another array. Each biological sample is hybridized exactly two times (Table 1.2).
 3.
Dyeswaps for which each couple of biological samples from two conditions are hybridized on two arrays with reverse dyes. Dyeswaps are constrained to be couplebalanced (Table 1.3).
Three different balanced reverse dye designs for the comparison of 2 treatments
1  array  1  2  3  4  5  6  7  8  9  10 

Cy5  A1  B5  A3  B9  A5  B6  A7  B10  A9  B9  
Cy3  B3  A2  B8  A4  B2  A6  B1  A8  B4  A10  
2  array  1  2  3  4  5  6  7  8  9  10 
Cy5  A1  B1  A2  B2  A3  B3  A4  B4  A5  B5  
Cy3  B1  A2  B2  A3  B3  A4  B4  A5  B5  A1  
3  array  1  2  3  4  5  6  7  8  9  10 
Cy5  A1  B1  A2  B2  A3  B3  A4  B4  A5  B5  
Cy3  B1  A1  B2  A2  B3  A3  B4  A4  B5  A5 
Dyeswap design is mostly used when the technical error is higher than the biological variability, either to reduce the technical variance, or when genespecific dyebias 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 Dyeswitch is used for the first and sometimes also for the second classes. Dyeswitch 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 Ttests (or regularized Ttests) 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 individuallybalanced dyeswitch design: two conditions A and B are compared in a twocolor 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 Cy 5 and the second labelled with Cy 3. Then 2n microarrays are hybridized with respectively A_{1}Cy 5 and B_{1}Cy 3, B_{1}Cy 5 and A_{2}Cy 3, A_{2}Cy 5 and B_{2}Cy 3, and so on till B_{ n }Cy 5 and A_{1}Cy 3, (see Table 1.2). There are 2n samples, 4n labelled samples, 2n microarrays, and each sample is hybridized two times (one with Cy 5 and one with Cy 3) 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 logscale, 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 genespecific. Here we use an analysis of variance (ANOVA) model for the expression measure as introduced by [10].
where

μ_{ A }and μ_{ B }are the population mean expression measures for condition A and B.

δ_{l(i)}is a twolevel 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 genespecific 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 Cy 3 and Cy 5 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 ${{T}^{\prime}}_{i}$ are independent random variables. T_{ i }and M_{ i }are the two components of the socalled technical error.
Model on the difference of expression on one array
which may be writtenD_{ i }= μ + BD_{ i }+ (1)^{i+1}δ + TD_{ i }
where

μ = μ_{ A } μ_{ B }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 $\sqrt{2}$σ_{ B },

TD_{ i }= T_{ i } ${{T}^{\prime}}_{i}$ is an independent random variable with mean 0 and standard deviation $\sqrt{2}$σ_{ T },

δ = δ_{1}  δ_{2} is the difference between the Cy3 and Cy5 dye effects. This term accounts for the genespecific dye bias.
with the convention that 2n + 1 = 1.
In this study, we present and compare different strategies for the statistical analysis of individuallybalanced 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 individuallybalanced designs. The details of the computation are given in the Methods section.
Results
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${T}_{N}=\sqrt{2}n\frac{\overline{D}}{\sqrt{{S}^{2}}}$
is computed.

Unbiased Paired Method (UP): for each gene, the unbiased paired statistic${T}_{UP}=\sqrt{2}n\frac{\overline{D}}{\sqrt{({S}^{2}+2C)}}$
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${T}_{UU}=\sqrt{n}\frac{\overline{X}\overline{Y}}{\sqrt{{S}_{X}^{2}+{S}_{Y}^{2}2{C}_{XY}}}$
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 ${S}_{X}^{2}+{S}_{Y}^{2}2{C}_{XY}$. Since C_{ XY }is nonnegative, 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).

Mixed Model with ML estimation (ML): for each gene, model (1) is adjusted with the Maximum Likelihood algorithm.

Mixed Model with REML estimation (REML): for each gene, model (1) is adjusted with the Restricted Maximum Likelihood algorithm.
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). We considered 3 different values for ${\sigma}_{B}^{2}$ (0.5, 1, 2) and ${\sigma}_{M}^{2}$ (1, 2, 5), 4 values for the number of samples in one condition (5, 10, 20, 30) and 5 possible values for the differential expression μ = μ_{ A } μ_{ B }(0, 1, 2, 3, 4). The parameter σ_{ T }was fixed at 1. For each combination of the parameters, 10,000 genes were simulated.
Control of the Type I error rate
Actual level of the 5 test procedures in one simulation of 10 000 genes
${\sigma}_{B}^{2}$ = 0.5  ${\sigma}_{B}^{2}$ = 2  

Method  5  10  20  30  5  10  20  30 
Naive  6.9 (0.2)  7.3 (0.2)  7.3 (0.2)  7.5 (0.2)  13.2 (0.3)  13.9 (0.3)  14.0 (0.3)  14.2 (0.3) 
Unbiased Paired  5.2 (0.2)  5.2 (0.2)  5.2 (0.2)  5.3 (0.2)  8.2 (0.3)  6.9 (0.2)  6 (0.2)  5.8 (0.2) 
Unbiased Unpaired  2.1 (0.1)  1.3 (0.1)  1.0 (0.1)  1 (0.1)  4.6 (0.2)  3.4 (0.2)  2.7 (0.1)  2.9 (0.2) 
ML  8.5 (0.3)  8.6 (0.3)  8.3 (0.3)  8.3 (0.3)  12.5 (0.4)  11.1 (0.3)  9.9 (0.3)  9.8 (0.3) 
REML  4.7 (0.2)  4.2 (0.2)  4.5 (0.2)  4.9 (0.2)  14.7 (0.4)  8.5 (0.3)  5.9 (0.2)  5.5 (0.2) 
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
Power of the UU, UP and REML test procedures
μ = 1  μ = 3  

Nb Samples  ${\sigma}_{B}^{2}$  UU  UP  REML  UU  UP  REML 
5  0.5  5.6  13.6  10.6  55.5  92.1  86.75 
5  2  2.8  5.0  12.95  17.9  29.4  34.75 
10  0.5  13.2  39.3  33.97  77.8  100.0  99.64 
10  2  3.5  7.8  9.06  45.0  63.5  63.06 
20  0.5  35.0  80.1  78.13  98.8  100.0  100.0 
20  2  7.3  14.5  13.93  82.6  94.8  94.53 
30  0.5  51.9  95.5  95.05  100.0  100.0  100.0 
30  2  12.1  22.5  21.74  96.2  99.6  99.53 
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 Ttest 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
CPU times of procedures UP and REML
n  UP CPU  REML CPU  No REML CV 

5  2.3  787  56.9 
10  2.6  212  5 
20  2.8  467  0 
30  3.2  1046  0.16 
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 dyeswitch 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.
Teleost fish dataset

N 1 → S 1 → N 3 → S 3 → N 5 → S 5 → N 2 → S 2 → N 4 → S 4 → N 1

S 1 → G 1 → S 3 → G 3 → S 5 → G 5 → S 2 → G 2 → S 4 → G 4 → S 1

N 1 → G 2 → N 3 → G 4 → N 5 → G 1 → N 2 → G 3 → N 4 → G 5 → N 1
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):Y_{ ijkg }= m + A_{ i }+ D_{ j }+ (AD)_{ ij }+ G_{ g }+ (AG)_{ ig }+ (DG)_{ jg }+ (V G)_{ kg }+ e_{ ijkg },
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.
Lists of genes for the Teleofish experiment
Oleksiak list [8]  UP list 

RAN GTP binding protein hypo P FLJ20727 ribosomal protein L27 dihydrolipoamide dehydrogenase GTP binding protein  
Steroidogenic acute regulatory protein hypo P FLJ11275 capping protein muscle Z line orla C4 surface glycoprotein HT7 precursor methionine adeno. regulatory Von Willebrand factor succinate dehydrogenase complex KIAA1481 protein protein disulfide isomerase annexin V  Thioredoxin nascent polypeptide associated dnaK type molec. chap. prec. ribosomal protein S16 
{S 1, S 2, S 4, S 5, N 1}, {G 1, G 2, G 3, G 4, G 5, S 3}, {N 2, N 5, N 3, N 4}.
These 3 classes roughly correspond to the three populations of interest, up to 2 misclassified observations. In the original article, the partition in 3 classes gave{N 1, N 2, N 3, N 4, N 5}, {S 1, S 4}, {G 1, G 2, G 3, G 4, G 5, S 2, S 3, S 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 dyeswitch experiments. We showed on simulations that the naive paired Ttest, 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.
For small sample size experiments, it is advised to use regularized Ttest, see [16–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 dyeswitch 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.
Methods
Paired test procedure
To perform a statistical test on parameter μ, we need to estimate ${V}_{\overline{D}}$.
Naive variance estimate
Unbiased variance estimate
which is approximately distributed as a Student distribution with 2n  2 df under H_{0}.
Unpaired test procedure
Declarations
Acknowledgements
The authors thank Douglas L. Crawford who provided the Teleofish dataset.
Authors’ Affiliations
References
 Yang Y, Speed T: Design issues for cDNA microarray experiments. Nat Rev Genet 2002, 3(8):579–88.PubMedGoogle Scholar
 Churchill G: Fundamentals of experimental design for cDNA microarrays. Nat Genet 2002, 32: 490–495. 10.1038/ng1031View ArticlePubMedGoogle Scholar
 Yang Y, Dudoit S, Luu P, Speed T: Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nuclear Acids Res 2002, 30(4):e15. 10.1093/nar/30.4.e15View ArticleGoogle Scholar
 Kerr M, Afshari C, Bennett L, Bushel P, Martinez J, Walker N, Churchill G: Statistical Analysis of a gene expression microarray experiment with replication. Statistica Sinica 2002, 12: 203–217.Google Scholar
 Dobbin K, Shih J, Simon R: Statistical design of reverse dye microarrays. Bioinformatics 2003, 19(7):803–810. 10.1093/bioinformatics/btg076View ArticlePubMedGoogle Scholar
 MartinMagniette M, Aubert J, Cabannes E, Daudin J: Evaluation of the genespecific dye bias in cDNA microarray experiments. Bioinformatics 2005, 21(9):1995–2000. 10.1093/bioinformatics/bti302View ArticlePubMedGoogle Scholar
 Dobbin K, Kawasaki E, Petersen D, Simon R: Characterizing dye bias in microarray experiments. Bioinformatics 2005, 21(10):2430–2437. 10.1093/bioinformatics/bti378View ArticlePubMedGoogle Scholar
 Oleksiak M, Churchill G, Crawford D: Variation in gene expression within and among natural populations. Nat Genet 2002, 32: 261–266. 10.1038/ng983View ArticlePubMedGoogle Scholar
 Whitehead A, Crawford D: Neutral and adaptive variation in gene expression. Proc Natl Acad Sci USA 2006, 103(14):5425–5430. 10.1073/pnas.0507648103PubMed CentralView ArticlePubMedGoogle Scholar
 Kerr M, Martin M, Churchill G: Analysis of variance for gene expression microarray data. J Comput Biol 2000, 7: 819–837. 10.1089/10665270050514954View ArticlePubMedGoogle Scholar
 Wit E, McClure J: Statistics for Microarrays: Design, Analysis and Inference. Chichester: Wiley; 2004.View ArticleGoogle Scholar
 Landgrebe J, Bretz F, Brunner E: Efficient twosample designs for microarray experiments with biological replications. Silico Biology 2004., 4:Google Scholar
 Landgrebe J, Bretz F, Brunner E: Efficient design and analysis of two colour factorial microarray experiments. Comput Stat & Data Anal 2006, 50(2):499–517. 10.1016/j.csda.2004.08.014View ArticleGoogle Scholar
 Mansouri N, Sandra O, Aubert J, Everts R, Galio L, Heyman Y, Audouart C, Degrelle S, Hue I, Yang X, Lewin H, JP R: Identification of differentially regulated genes in the endometrium of cyclic and pregnant cows using a highthroughput transcriptome analysis. American Journal of Reproductive Immunology 2007, 58(5):204–220.Google Scholar
 Wit E, Nobile A, Khanin R: Nearoptimal designs for dual channel microarray studies. J R Stat Soc C 2005, 54(5):817–830. 10.1111/j.14679876.2005.00519.xView ArticleGoogle Scholar
 Baldi P, Long A: A Bayesian Framework for the Analysis of Microarray Expression Data: Regularized tTest and Statistical Inferences of Gene Changes. Bioinformatics 2001, 17: 509–519. 10.1093/bioinformatics/17.6.509View ArticlePubMedGoogle Scholar
 Delmar P, Robin S, Daudin J: VarMixt: efficient variance modelling for the differential analysis of replicated gene expression data. Bioinformatics 2005, 21(4):502–508. 10.1093/bioinformatics/bti023View ArticlePubMedGoogle Scholar
 Smyth G: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 2004, 3: Article 3. 10.2202/15446115.1027View ArticleGoogle Scholar
 Tusher V, Tibshirani R, Chu C: Significance analysis of microarrays applied to transcriptional response to ionizing radiations. PNAS 2001, 98: 5116–5121. 10.1073/pnas.091062498PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.