A power law global error model for the identification of differentially expressed genes in microarray data
 Norman Pavelka†^{1},
 Mattia Pelizzola†^{1},
 Caterina Vizzardelli†^{1},
 Monica Capozzoli^{1},
 Andrea Splendiani^{1},
 Francesca Granucci^{1} and
 Paola RicciardiCastagnoli^{1}Email author
DOI: 10.1186/147121055203
© Pavelka et al; licensee BioMed Central Ltd. 2004
Received: 16 July 2004
Accepted: 17 December 2004
Published: 17 December 2004
Abstract
Background
Highdensity oligonucleotide microarray technology enables the discovery of genes that are transcriptionally modulated in different biological samples due to physiology, disease or intervention. Methods for the identification of these socalled "differentially expressed genes" (DEG) would largely benefit from a deeper knowledge of the intrinsic measurement variability. Though it is clear that variance of repeated measures is highly dependent on the average expression level of a given gene, there is still a lack of consensus on how signal reproducibility is linked to signal intensity. The aim of this study was to empirically model the variance versus mean dependence in microarray data to improve the performance of existing methods for identifying DEG.
Results
In the present work we used data generated by our lab as well as publicly available data sets to show that dispersion of repeated measures depends on location of the measures themselves following a power law. This enables us to construct a power law global error model (PLGEM) that is applicable to various Affymetrix GeneChip data sets. A new DEG identification method is therefore proposed, consisting of a statistic designed to make explicit use of modelderived measurement spread estimates and a resamplingbased hypothesis testing algorithm.
Conclusions
The new method provides a control of the false positive rate, a good sensitivity vs. specificity tradeoff and consistent results with varying number of replicates and even using single samples.
Background
DNA microarrays have become common tools for monitoring genomewide expression in biological samples harvested under different physiological, pathological or pharmacological conditions. One of the most challenging problems in microarray data analysis is probably the identification of differentially expressed genes (DEG) when comparing distinct experimental conditions. In spite of its biological relevance, there is still no commonly accepted way to answer this question.
An ideal DEG identification method should limit both false positives, i.e. genes wrongly called significant (type 1 errors), and false negatives, i.e. genes wrongly called not significant (type 2 errors). To this end, understanding how gene expression values measured in replicated experiments are spread around the true expression level of each gene, would help to distinguish biologically relevant gene expression changes from fluctuations due to different sources of variability that are unrelated to the biological phenomenon under investigation. Measurement error estimates can be obtained in two ways: either by empirically inferring noise from highly replicated data or by deducing noise from a theoretical error model [1]. Especially when the experimental design requires the investigation of a high number of conditions, the former strategy is not always feasible, because of the high cost of these experiments or due to the availability of biological material. In addition, there is still a lack of consensus on how gene expression values from replicated experiments should be theoretically distributed, which restricts the application also of the latter strategy.
The most widely used methods for identifying DEG range from purely empirical filtering techniques (e.g. selecting genes that show a fold change higher than a fixed threshold) to more sophisticated statistical tests such as the signaltonoise ratio described by Golub et al. [2] or the Significance Analysis of Microarrays (SAM) method by Tusher et al. [3]. While empirical filtering techniques rely on arbitrarily chosen thresholds and are unable to provide any type of control on the significance of the results, the more sophisticated statistical tests usually need a high degree of replication in the data to accurately measure genespecific variability.
In the past years various authors have proposed competing error models for microarray data from which discordant implications for the variance versus mean dependence can be deduced. Chen et al. [4] first proposed a simple Gaussian model, more recently Ideker et al. [5] and Li and Wong [6] introduced twocomponent models containing a multiplicative and an additive error term. All of these models implicitly or explicitly assume a constant coefficient of variation (CV), implying that standard deviation should vary proportionally with the mean. More recently, Rocke and Durbin [7] proposed a variation of the twocomponent model from which they derived that variance of repeated microarray measures is a quadratic function of the mean. Dealing specifically with spotted cDNA microarray technology Baggerly et al. [1] proposed a betabinomial model, from which it can be derived that variance is a secondorder polynomial function of the mean. Unfortunately, most of these models are based on theoretical assumptions that have been verified on simulated data or on data sets consisting of small numbers of replicates. More recently, Tu et al. [8] empirically modeled the variance versus mean dependence from a data set consisting of ten replicated oligonucleotide microarray experiments. According to the authors, the variance of the genes should decay exponentially with the mean, but only for moderately expression values. Taken together, all these aspects could limit the applicability of these error models.
Independently from the choice of the error model, another point that remains to be faced is on how to estimate residual error. A discussed by Wright et al. [9] the possibilities range between two extremes: either obtaining a single variance estimate across all genes or obtaining a genespecific residual variance. In the same paper a hybrid approach is proposed in which information from all genes is used to fit a single linear model from which the genespecific variance estimates can be deduced. In the present work we chose to follow an approach similar to the latter.
The aim of this study was to use highly replicated microarray data to empirically determine the true variance versus mean dependence that exists in this type of data. This knowledge enabled the proposal of PLGEM as a simple but powerful error model. We fitted the proposed model on various data sets without prefiltering the data, deriving an improved test statistics and identifying DEG even in data sets with very little number of replicates.
Results
Variance versus mean dependence
The power law global error model
Based on the previous observation, we chose to empirically model measurement noise through linear regressions:
where s and respectively represent standard deviation and mean of repeated measures. Error term ε is the realization of a random variable E that we will show later to be normally distributed as assumed when fitting a linear model. Inspired by the previous experimental observations we propose the following power law global error model (PLGEM):
Model parameters α and β can be estimated from linear regression coefficients in 1 in a straightforward way:
α = e^{ c } eqn. 3
β = k eqn. 4
PLGEM fitting method
Instead of performing a simple linear fit through the whole set of points, we preferred to implement a method that could provide improved model robustness by partitioning the data to gain local estimates of spread as in Mutch et al. [9]. Most importantly, this method should also provide the possibility to choose different levels of confidence when modeling the spread of the data. Note that Mutch et al. [9] proposed to model withinreplicates fold changes as a function of average expression using a model that was very different from PLGEM. Therefore, the following algorithm was applied:

Finally, find a linear fit through the set of p modeling points using the leastsquares method and obtain a slope k and an intercept c of the resulting regression function.
Effect of the presence of DEG when applying the classic permutation strategy to the PLGEMSTN statistic.
FPR vs. significance level  estimated vs. observed FDR  

including DEG  excluding DEG  including DEG  excluding DEG  
# of genes in data set  slope  intercept  adj. R^2  slope  intercept  adj. R^2  slope  intercept  adj. R^2  slope  intercept  adj. R^2 
22300  1.690  3.070  0.856  0.871  0.114  0.938  0.187  0.383  0.314  1.090  0.692  0.826 
10000  1.710  2.470  0.881  0.888  0.083  0.931  0.224  0.247  0.433  1.040  0.555  0.824 
5000  1.460  1.270  0.880  0.877  0.147  0.934  0.093  0.228  0.067  1.080  0.425  0.836 
2500  1.590  1.180  0.888  0.864  0.155  0.939  0.092  0.176  0.135  1.110  0.348  0.853 
1500  1.670  0.935  0.908  0.876  0.166  0.948  0.078  0.122  0.170  1.130  0.182  0.880 
1000  1.720  0.689  0.874  0.944  0.002  0.956  0.038  0.122  0.082  0.991  0.226  0.897 
500  1.900  0.263  0.864  0.857  0.255  0.956  0.062  0.030  0.307  1.160  0.038  0.909 
200  2.490  0.059  0.875  0.905  0.378  0.946  0.064  0.012  0.426  1.050  0.233  0.914 
As a straightforward application of this modeling method, PLGEM could be fitted at the 50^{th}percentile to obtain a central tendency of standard deviation to be used for improving test statistics (see next section). Another application of this method could be to fit PLGEM at the 5^{th} and at the 95^{th}percentile of standard deviation, to consequently find the limits of the corresponding 90% empirical confidence interval of standard deviation.
In order to verify the feasibility of the former application, fitting of PLGEM on reallife data as well as distribution properties of the random variable E were investigated by analyzing the residuals of the model, i.e. differences between observed and expected values:
Improved teststatistics for detecting differential expression
In order to identify DEG, we implemented the following general algorithm derived from the framework of statistical hypothesis testing, in which we test against the null hypothesis of nondifferential expression. First of all, we chose to implement as the test statistic the signaltonoise ratio (STN) already used by Golub et al. [2], because it explicitly takes unequal variances into account and because it penalizes genes that have higher variance in each class more than those genes that have a high variance in one class and a low variance in another [11]:
Identification of differentially expressed genes
A resamplingbased method for estimating the null distribution
Though ranking of genes based on the absolute value of their teststatistic has been proven to be an effective method for selecting DEG, an even more useful way would be to compare the observed statistic with its null distribution (the distribution of values of the statistic that are expected by chance for a not differentially expressed gene), in order to control the FPR.
A classic approach to empirically obtain the null distribution of a teststatistic is running a series of random permutations of the chip indexes of the full data set and recomputing the teststatistics at each permutation. Permutated teststatistics can then be pooled and significance thresholds (i.e. expected false positive rates) are found as specific quantiles of the null distribution.
Nevertheless, we can foresee that the classic permutation strategy may not be optimal for estimating the actual FPR when the teststatistic makes use of a global error model such as PLGEM. We can in fact hypothesize that measurement spread of DEG may not be accurately described by means of a global error model that was designed to describe signal variability in absence of differential expression. To test this hypothesis we compared the correlation between the expected significance level and the observed FPR using PLGEMSTN and the classic permutation strategy either including or excluding DEG during the permutation step. To this end, data sets containing different percentages of DEG were obtained by merging the 62 known DEG of the Latin Square data set with differently sized random samples of not DEG extracted from the same data set. As predicted, the presence of DEG during the permutation step caused the significance level to be less correlated with the observed FPR and this correlation worsened with increasing percentages of DEG (Table 1). This lack of correlation was dramatically amplified when expected and observed numbers of false positives were divided by the number of selected genes to obtain an oversimplified estimate of the false discovery rate (FDR) and the observed FDR. Conversely, when DEG were omitted during the permutation step the correlation between estimated and observed FPR or estimated and observed FDR was sensibly higher for each tested percentage of DEG. We hereby by no means claim that this FDR estimate is the most accurate. A more appropriate relationship between FPR and FDR can be found in the paper by Storey and Tibshirani [13]. Nevertheless, the explicit control of the FDR goes beyond the scope of the present paper.
Since in reallife data sets true DEG are unknown in advance, we propose the following resamplingbased method to obtain the null distribution of not DEG when comparing n_{1} replicates of condition A with n_{2} replicates of condition B:

Artificial condition A* is obtained by randomly sampling with replacement n_{1} indexes corresponding to the replicates of only one experimental condition. If available, chose the condition with the highest number of replicates;

Similarly sample n_{2} values from the same set to obtain indexes of artificial condition B*;

Compute resampled teststatistics between A* and B* at each cycle.
In accordance with the previous observations, the ROC curve of the resampling method applied to the PLGEMSTN statistic was not significantly different from the ROC curve of the classic permutation strategy (excluding DEG) applied to the same statistic on the Latin Square data set (data not shown). Conversely, ROC curves of the classic permutation strategy (including DEG) applied to the CLASSICSTN statistic and of the SAM method gave poorer performance similarly to the results in Figure 3 (data not shown).
Increased robustness to varying number of replicates
We finally evaluated the possibility of applying our method also to data sets where one of the experimental conditions was investigated only with a single sample without replication. To this end, we used the remaining four reduced data sets that could not be used in the previous comparison. In this case, the same PLGEM parameters derived from the sixteen baseline columns were applied to each of the single LPStreated DC sample to obtain an estimate of standard deviation associated to each gene expression value, treated here as if it was a mean value from a larger group of values. Interestingly, when results obtained through this procedure were compared to the previously described results a comparable number of DEG was identified and only one probe set was newly detected in comparison to the previously identified ones (data not shown), arguing for a good consistency of results.
Discussion
PLGEM accurately describes GeneChip data variability
In the present work we described a new global error model for microarray gene expression data that describes measurement variability with the same degree of accuracy over the whole dynamic range of values and that can be fitted at any desired quantile of spread. PLGEM has proven to correctly model signal standard deviation, in spite of the presence of different sources of variability, e.g. biological variability as well as the use of different target preparation protocols or of different chips. Moreover, PLGEM has shown to be able to deal with the great variability that exists at low expression levels while at the same time considering the significant relative reproducibility of highly expressed genes. Previously proposed error models assumed that measurement spread depended on signal location following different mathematical relationships, but none of them was based on a power law thus far. Analysis of the residuals showed a good fit of PLGEM to a number of highdensity oligonucleotide microarray data sets, with model parameters being very similar to each other even when dealing with RNA samples coming from completely different biological sources and analyzed on different array layouts. This suggests that PLGEM could represent a general Affymetrix GeneChip measurement noise model. Even though scaled MAS5 Signals gave satisfactory modeling results, a further improvement could be achieved by using other emerging gene expression indices [6, 14] or more sophisticated normalization techniques, e.g. quantile normalization [15]. Interestingly, if the same evaluation of sensitivity vs. specificity using ROC plots on the Latin Square data set was done using GCRMA expression values [16], the results were even more striking than using MAS5 Signals (data not shown). Further studies will be needed to assess if PLGEM is also able to deal with data coming from microarray technologies others than Affymetrix GeneChips.
Interestingly, model parameter β was found to be quite stable and comprised between 0 and 1 in all analyzed data sets. It is noteworthy that for β ∈ (0:1) absolute variability increases with growing expression values, while relative variability decreases (compare panel B with panel D of Figure 1). On the other hand, none of the models mentioned in the background section seem to agree with these experimental observations. Formal statistical reasoning could unravel the underlying theoretical error model that leads to the power law relationship that was observed to be at the basis of the variance versus mean dependence in replicated microarray data.
A PLGEMbased method successfully detects differential expression
In spite of the lack of a theoretical statistical model, the empirical model presented here has proven its applicability in the identification of DEG, providing improved results under a wide range of different testing conditions. In comparison to other commonly used DEG identification methods, the proposed approach demonstrated improved specificity and sensitivity on the Latin Square data set and robustness to decreasing number of replicates on the 16iDC+LPS data set. The good performance of our proposed method is reasonably due to the fact that it relies on a global error model. As an example, when the classic permutation strategy is applied to the CLASSICSTN statistic or when the SAM method is used, the selected genes are apparently more dependent on the number and identity of the replicates than when our proposed approach is used. We hypothesize that, when no error model is assumed and a small number of replicates is present in the data set, the probability of observing for some genes coincidently very similar (or very dissimilar) values increases, thus leading to an underestimation (or overestimation) of the standard deviation and a consequent overestimation (or underestimation) of the test statistic, finally leading to false positives (or false negatives).
Interestingly, when the performance of our method was compared on a data set of DC stimulated for 24 hours with LPS, SAM showed a decreased sensitivity in identifying downregulated genes when the number of LPS replicates was low (data not shown). Under these experimental conditions DC undergo a process known as maturation, which is a specialized form of cellular differentiation, for which both up and downregulation of gene expression is expected [17, 18]. We speculate that SAM did not select these genes, because of the combination of two effects. First of all, downregulated genes are expected to have lower and therefore intrinsically more variable expression values in the four LPS replicates than in the sixteen replicates of immature DC. When, in addition, the number of LPS replicates becomes too low, SAM filters these genes out to control the FDR. In agreement with this hypothesis SAM was perfectly able to identify downregulation when the full data set was used (data not shown).
The gene selection method proposed in the present work does not provide a direct control on the FDR, but the significance level has been proven to be a direct estimate of the FPR. Thus, if a significance level of 0.001 is used and 12488 probe sets are displayed on the MGU74Av2 chip, 12–13 genes are expected to be selected by chance in cases where all genes are in fact not differentially expressed. Therefore, a researcher can test how many genes would be selected over a range of different significance levels and chose the one that results in the most acceptable compromise between number of selected genes and estimated FPR.
Conclusions
The proposed DEG identification method provides a direct control of the FPR and an indirect control of the FDR. Moreover, as tested on the Latin Square data set, our method improved the specificity vs. sensitivity tradeoff in comparison to other commonly applied DEG selection techniques. It finally showed an increased robustness when different replicates or numbers thereof are analyzed, giving consistent results even in data sets containing single samples. In conclusion, the global error model presented here may facilitate the analysis of microarray gene expression data by discriminating information from noise, and thus possibly helping the formulation of new hypothesis concerning gene functions.
Methods
Data sets
16iDC
RNA was harvested from ten biological samples of unstimulated immature mouse dendritic cells (DC), each extracted from an independent batch of cells. One operator prepared the biotinlabeled cRNA for hybridization from three of the ten RNA samples, a second operator prepared the remaining seven. While operator 1 applied the total RNA protocol to all of its three samples, operator 2 applied the purified mRNA protocol to five of its seven samples and the total RNA protocol to the remaining two. Two of the three cRNA samples prepared by operator 1 and four of the seven cRNA samples prepared by operator 2 have been hybridized twice; therefore, a total of 16 MGU74Av2 GeneChips (Affymetrix, Santa Clara, CA) have been employed.
Leigh syndrome
Eight RNA samples were harvested from human fibroblast cell lines each deriving from a distinct Leigh syndrome patient [19, 20] and individually hybridized on HGU133A GeneChips (Affymetrix).
Muscle biopsies
Four individual and two pooled RNA samples from human muscle biopsies of sixteen healthy young male donors were hybridized on six HGU95Av2 GeneChips (Affymetrix). This data set was downloaded from [21], experiment code: GSE80 [22].
Latin Square
This data set consists of 3 technical replicates of 14 separate hybridizations (named Exp01–14) of 42 spiked transcripts in a complex human background at concentrations ranging from 0.125 pM to 512 pM. Thirty of the spikes are isolated from a human cell line, four spikes are bacterial controls, and eight spikes are artificially engineered sequences believed to be unique in the human genome. Further details on the design of the Latin Square data set can be found at [23]. Considering the redundancy of some probe sets, there are a total of 62 distinct probe sets designed to match the 42 spiked transcripts.
16iDC+LPS
This data set consists of the same samples of the 16iDC data set, but includes additional four samples as a second experimental condition. To this end dendritic cells were stimulated to mature with lipopolysaccharide (LPS) for 24 hours. Two independent biological samples were harvested and individually processed by the same two operators that prepared the samples for the 16iDC data set: one applied the total RNA protocol, the other one applied the purified mRNA protocol. Each cRNA sample was hybridized twice, thus using a total of four Affymetrix MGU74Av2 chips.
Software
All chips mentioned in the present study were hybridized and scanned following Affymetrix recommendations and MicroArray Suite 5.0 (MAS5) was used as the image acquisition and analysis software. All data sets used passed quality control tests and probe set signals were scaled so that the 4%trimmed mean of all expression values of each chip was equal to a predefined reference intensity (called TGT) following manufacturer's recommendations:
TGT = 100 for MGU74Av2 and HGU133A chips and TGT = 500 for HGU95Av2 chips.
All procedures for fitting PLGEM, for calculating observed PLGEMbased signaltonoise ratios (STN), for obtaining expected PLGEMSTN through the resamplingbased approach and for comparing observed with expected STN values have been implemented as R functions [24] and will be soon submitted for integration into the Bioconductor project [25].
Notes
Declarations
Acknowledgements
This research was supported by AIRC (Italian Association for Cancer Research) and by a grant from the CARIPLO Foundation ("Development of Functional Genomics and Bioinformatics platforms aimed to foster novel approaches in immunotherapy and molecolar diagnostics"). Authors would like to kindly acknowledge Stefano Monti (Broad Institute, Cambridge, MA, USA) and Peter J. Park (Harvard Medical School, Boston, MA, USA) for helpful discussion and Ottavio Beretta, Gianpiero Cattaneo, Maria Foti, Giorgio Moro and Ettore Virzi (University of MilanoBicocca, Milan, Italy) for critical reading of the manuscript. We are also grateful to Valeria Tiranti, Rossana Mineri and Massimo Zeviani (Unit of Molecular Neurogenetics, National Neurological Institute "Carlo Besta", Milano, Italy) for kindly providing the HGU133A Leigh syndrome data set.
Authors’ Affiliations
References
 Baggerly KA, Coombes KR, Hess KR, Stivers DN, Abruzzo LV, Zhang W: Identifying differentially expressed genes in cDNA microarray experiments. J Comput Biol 2001, 8: 639–659. 10.1089/106652701753307539View ArticlePubMedGoogle Scholar
 Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA, Bloomfield CD, Lander ES: Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science 1999, 286: 531–537. 10.1126/science.286.5439.531View ArticlePubMedGoogle Scholar
 Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A 2001, 98: 5116–5121. 10.1073/pnas.091062498PubMed CentralView ArticlePubMedGoogle Scholar
 Chen Y, Dougherty ER, Bittner ML: Ratio based decisions and the quantitative analysis of cDNA microarray images. J Biomed Opt 1997, 2: 364–374. 10.1117/1.429838View ArticlePubMedGoogle Scholar
 Ideker T, Thorsson V, Siegel AF, Hood LE: Testing for differentiallyexpressed genes by maximumlikelihood analysis of microarray data. J Comput Biol 2000, 7: 805–817. 10.1089/10665270050514945View ArticlePubMedGoogle Scholar
 Li C, Wong WH: Modelbased analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Natl Acad Sci U S A 2001, 98: 31–36. 10.1073/pnas.011404098PubMed CentralView ArticlePubMedGoogle Scholar
 Rocke DM, Durbin B: A model for measurement error for gene expression arrays. J Comput Biol 2001, 8: 557–569. 10.1089/106652701753307485View ArticlePubMedGoogle Scholar
 Tu Y, Stolovitzky G, Klein U: Quantitative noise analysis for gene expression microarray experiments. Proc Natl Acad Sci U S A 2002, 99: 14031–14036. 10.1073/pnas.222164199PubMed CentralView ArticlePubMedGoogle Scholar
 Wright GW, Simon RM: A random variance model for detection of differential gene expression in small microarray experiments. Bioinformatics 2003, 19: 2448–2455. 10.1093/bioinformatics/btg345View ArticlePubMedGoogle Scholar
 Mutch DM, Berger A, Mansourian R, Rytz A, Roberts MA: The limit fold change model: A practical approach for selecting differentially expressed genes from microarray data. BMC Bioinformatics 2002, 3: 17. 10.1186/14712105317PubMed CentralView ArticlePubMedGoogle Scholar
 Cancer Genomics Software: GeneCluster 2.0b Reference Manual[http://www.broad.mit.edu/cancer/software/genecluster2/gc_ref.html]
 Broberg P: Statistical methods for ranking differentially expressed genes. Genome Biol 2003., 4: epub May 29Google Scholar
 Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci U S A 2003, 100: 9440–9445. 10.1073/pnas.1530509100PubMed CentralView ArticlePubMedGoogle Scholar
 Irizarry RA, Hobbs B, Collin F, BeazerBarclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 2003, 4: 249–264. 10.1093/biostatistics/4.2.249View ArticlePubMedGoogle Scholar
 Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on bias and variance. Bioinformatics 2003, 19: 185–193. 10.1093/bioinformatics/19.2.185View ArticlePubMedGoogle Scholar
 Wu Z, Irizarry RA, Gentleman R, Murillo FM, Spencer F: A Model Based Background Adjustment for Oligonucleotide Expression Arrays. Johns Hopkins University, Dept. of Biostatistics Working Papers, Working Paper 1 (May 28, 2004) [http://www.bepress.com/jhubiostat/paper1]
 Granucci F, Vizzardelli C, Virzi E, Rescigno M, RicciardiCastagnoli P: Transcriptional reprogramming of dendritic cells by differentiation stimuli. Eur J Immunol 2001, 31: 2539–2546. 10.1002/15214141(200109)31:9<2539::AIDIMMU2539>3.3.CO;20View ArticlePubMedGoogle Scholar
 Granucci F, Vizzardelli C, Pavelka N, Feau S, Persico M, Virzi E, Rescigno M, Moro G, RicciardiCastagnoli P: Inducible IL2 production by dendritic cells revealed by global gene expression analysis. Nat Immunol 2001, 2: 882–888. 10.1038/ni0901882View ArticlePubMedGoogle Scholar
 Tiranti V, Hoertnagel K, Carrozzo R, Galimberti C, Munaro M, Granatiero M, Zelante L, Gasparini P, Marzella R, Rocchi M, BayonaBafaluy MP, Enriquez JA, Uziel G, Bertini E, DionisiVici C, Franco B, Meitinger T, Zeviani M: Mutations of SURF1 in Leigh disease associated with cytochrome c oxidase deficiency. Am J Hum Genet 1998, 63: 1609–1621. 10.1086/302150PubMed CentralView ArticlePubMedGoogle Scholar
 Tiranti V, Jaksch M, Hofmann S, Galimberti C, Hoertnagel K, Lulli L, Freisinger P, Bindoff L, Gerbitz KD, Comi GP, Uziel G, Zeviani M, Meitinger T: Lossoffunction mutations of SURF1 are specifically associated with Leigh syndrome with cytochrome c oxidase deficiency. Ann Neurol 1999, 46: 161–166. Publisher Full Text 10.1002/15318249(199908)46:2<161::AIDANA4>3.0.CO;2OView ArticlePubMedGoogle Scholar
 Gene Expression Omnibus database[http://www.ncbi.nlm.nih.gov/geo]
 Welle S, Brooks AI, Thornton CA: Computational method for reducing variance with Affymetrix microarrays. BMC Bioinformatics 2002, 3: 23. 10.1186/14712105323PubMed CentralView ArticlePubMedGoogle Scholar
 Affymetrix – Latin Square Data[http://www.affymetrix.com/support/technical/sample_data/datasets.affx]
 The R Project for Statistical Computing[http://www.rproject.org]
 Bioconductor[http://www.bioconductor.org]
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.