Incorporation of gene-specific variability improves expression analysis using high-density DNA microarrays
© Budhraja et al; licensee BioMed Central Ltd. 2003
Received: 18 August 2003
Accepted: 28 November 2003
Published: 28 November 2003
The assessment of data reproducibility is essential for application of microarray technology to exploration of biological pathways and disease states. Technical variability in data analysis largely depends on signal intensity. Within that context, the reproducibility of individual probe sets has not been hitherto addressed.
We used an extraordinarily large replicate data set derived from human placental trophoblast to analyze probe-specific contribution to variability of gene expression. We found that signal variability, in addition to being signal-intensity dependant, is probe set-specific. Importantly, we developed a novel method to quantify the contribution of this probe set-specific variability. Furthermore, we devised a formula that incorporates a priori-computed, replicate-based information on probe set- and intensity-specific variability in determination of expression changes even without technical replicates.
The strategy of incorporating probe set-specific variability is superior to analysis based on arbitrary fold-change thresholds. We recommend its incorporation to any computation of gene expression changes using high-density DNA microarrays. A Java application implementing our T-score is available at http://www.sadovsky.wustl.edu/tscore.html.
The introduction of microarray technology has enabled investigators to profile the expression of a large number of genes, derived from diverse biological conditions, in a single experiment. Nevertheless, these experiments are expensive, and the cost is amplified by replication when data reliability is lessened. Technical and biological variability, key determinants of microarray reliability, are critical for assessing which genes are differentially expressed. The level of expression of gene products is estimated using a set of oligonucleotide probes (termed here "probe set"). Within the context of technology-related variability, the reproducibility of individual probes has been insufficiently addressed. [1–6]. Importantly, the impact of probe-specific reliability on data reproducibility, which directly influences experimental design, data analysis and interpretation, remains largely unexplored. This problem is further amplified with the use of oligonucleotide microarrays, which may be more susceptible to probe-specific variability than spotted cDNA arrays [7–10].
Static fold change metrics as an unbiased predictor of differentially expressed genes depend on the assumption of constant coefficient of variation . Since violation of this assumption is relatively common in microarray data, many methods have been designed to circumvent its requirement [2, 12–14]. These methods model variance as a function of intensity, and assume independence of probe sets. Other methods, including popular permutation tests and t-tests that do not explicitly rely on this assumption require replication in order to assess the variability associated with each probe set [15–17]. We proposed to use intensity-corrected measures of variance and a correlation test to determine if the assumption of probe set and variance independence is valid. Using a large replicate data set we found that signal variability in microarray data is in fact probe set-specific, and developed a novel method to integrate a priori-generated, replicate-based information on signal intensity and probe set variability into profiling differential gene expression. The signal intensity as well as probe set information serves as a database for future experiments performed without technical replicates. The use of our method enhances the significance of differences found using reliable probe sets and diminishes the significance of differences found using unreliable probe sets.
Results and Discussion
We have developed a large replicate data set, based on gene expression in a single pool of primary placental trophoblast cells. The cells were divided to three groups and exposed to two different peroxisome proliferator activated receptor gamma (PPARγ) ligands, troglitazone or GW7845, or to control, as described in Methods. Prior to hybridization each labeled RNA sample was divided to five aliquots. Each of the five cRNA aliquots was sequentially hybridized to identical lot number U95A, U95B, U95C, U95D, and U95E arrays, resulting in expression data for approximately 60 000 genes for each of the three conditions, for a total of 180 000 probe sets, each sampled in five replicates using a total of 75 chips.
The five replicates for each experimental condition i and probe set j yielded a mean ( ) and standard deviation ( ). We formulated an estimate of the standard deviation ( ) as a function of signal intensity using locally weighted scatter smooth plot (LOESS) local regression (using the PROC LOESS command in SAS). LOESS regression restricts attention to a small window of the data and fits a regression to that data. The window is then shifted and another local regression is calculated. These sequential regressions are combined to yield a LOESS curve [18–20]. Thus the LOESS curve indicates the average / estimated standard deviation associated with any given mean intensity, based on the replicate data set of ~180 000 sets of five replicates. To find a probe-specific effect, we compared the observed standard deviation to the estimated standard deviation. We define α i,j as:
where the function f is the LOESS regression, which returns the average standard deviation ( ) for a given intensity level ( ). By using the ratio of observed standard deviation to estimated standard deviation, α i,j represents a measure of residual variance for each probe set after correction for intensity.
Correlation of residual variance values across all three paradigms1.
Correlation of Alpha
Standard Deviation (after transformation 2 )
Control / GW7845
P < 0.0001
P < 0.0001
P < 0.0001
Control / Troglitazone
P < 0.0001
P < 0.0001
P < 0.0001
GW7845 / Troglitazone
P < 0.0001
P < 0.0001
P < 0.0001
To incorporate intensity and probe set-specific variability into determination of differentially expressed genes, we introduce novel methodology and an ad hoc t-statistic, called T-score. We have previously shown  that for two conditions i = 1 and i = 2 on probe set j, our intensity corrected T-score is derived as follows:
Where f(z j ) gives the predicted standard deviation value corresponding to the observed mean intensity of x 1,j and x 2,j derived from the LOESS regression. This T-score was shown to be independent of signal intensity. To correct for probe set-specific variability we multiplied the predicted standard deviation by (since α i,j is a ratio, we define to be the geometric mean of α i,j across all i paradigms) and obtained an intensity- and gene-specific estimate for the underlying variance. This makes sense intuitively, since α i,j is the ratio of observed to expected standard deviation, and is the mean of α i,j over the three paradigms. We noted that each of the 60 000 calculations of is based on only 15 observations. Therefore, a few observations for each probe set may inevitably result in falsely deviant values for . Cognizant of this possibility, we diminished the effect of extreme outliers of by using in its place. This function has the convenient properties of being symmetric with respect to the geometric mean and bringing extreme outliers closer to the null effect of one. The combined T-score becomes:
The replicate-based a priori generation of creates a database of probe set-specific variability coefficients that can be used for computation of T-score even from data derived without replicates, as conducted by most researchers. Instead of estimating variance with costly replicates, we contend that similar results can be obtained using estimates of variance generated from the intensity and probe set database we have created. Nevertheless, this simplification is not necessary. If an experiment has replicates the sample sizes can be carried through for their effect on the distribution of the mean. The T-score becomes:
and n 1 and n 2 represent the number of observations in paradigm 1 and 2, respectively. Here, we would disregard the observed variances and use the intensity and probe reliability index to create an estimate of the variance. A Bayesian approach could also be used with arbitrary degrees of freedom to incorporate the observed sample variances.
Correlation of predictive statistics with gene expression changes, determined by real-time quantitative RT-PCR.
Method of Correlation
T no gene correction
Generation of replicate data in microarray experiments can be utilized to assess variability, and consequently enhance data consistency . Whereas our approach to identify bias related to intensity-dependent variance is consistent with that of others [1, 2, 13, 15, 23, 24], we have also shown that consideration of probe set-specific variance is critical, given that the often-used assumption of probe set independence is false. Practically, we demonstrated that small expression differences detected by a more reliable probe set might be more important than larger expression differences detected by a less reliable probe set. Therefore, our analysis indicates that an optimal method for determining differentially expressed genes must account and correct for reliability of each individual probe set. While previous efforts have used replication to distinguish non-functional, or highly unreliable probe sets for elimination [1–5], none has incorporated probe set-specific reliability into expression change statistics. Our methodology is the first that allows independent use of previously obtained information on probe set reliability in subsequent experiments. Therefore, inference from prior estimates of probe set-specific variances into new experiments could not have been utilized using previously published approaches.
We were reassured by the fact that the correlation of our T-score method was similar to that of the t-test based Cyber-T . However, while exhibiting a similar performance, Cyber-T does not incorporate an a priori-defined gene correction factor, and relies on replication for estimation of variance. In contrast, our T-score approach integrates previously defined probe set-specific variance (via our database of α j ), defined by means of additional experimental paradigms (e.g., other ligands) for the same probe sets. Thus T-score, which independently accounts for intensity and probe set variances, may be utilized in array experiments even when performed without replication. Nevertheless, when adequate replicates are available, t-test is a suitable approach. It is also important to note that our analysis exclusively focuses on technical variability. While we demonstrate that a previously defined probe-specific variance can substitute for technical replicates, biological replicates are paramount to enhance accuracy of microarray-based expression analysis. In addition, the p-values associated with our empirically derived T-score do not include a correction for multiple comparisons, which should be accounted for when comparisons of expression level among thousands of genes are made [25, 26].
Another fundamental strength of our approach is the novel analysis of an extraordinarily large data set. Although our downstream analysis demonstrated the superiority of our methodology over an arbitrary cut-off approach, our results are limited by the fact that a fraction of our analyzed gene pool was expressed at a low level across all three paradigms. This might have led to an erroneous estimation of the specific variability of some probe sets. It should also be noted that our analysis is based on Affymetrix U95A-E gene-chip microarrays. Affymetrix has recently generated a new chip set (U133). While this may represent a technological advancement, the principles underlying our novel methodology are not addressed by the new chip-set. Whereas T-score values reported here are applicable only to the U95A-E set, the basis and principles underlying our analysis are applicable to any oligonucleotide microarray. We not only provide the first definitive proof that probe set-specific variability exists, but also offer the first generic methodology designed to utilize this information without performing costly replicate experiments. Probe reliability information based on cRNA targets that are expressed at high and low levels could be generated by biotechnology companies specializing in microarrays. This information, provided in conjunction with commonly available changes in p-value and determination of transcript presence, may serve to correct for technical variability in array experiments performed without replication.
We have developed a large replicate data set, based on gene expression in a single pool of primary trophoblast cells. Procurement of the placentas used in this study was approved by the human studies committee at Washington University School of Medicine, St. Louis, Missouri, USA. Primary human trophoblasts were prepared from three normal term human placentas as previously described  with previously published modifications , and cultured in 10 cm plates as previously detailed by our lab . Four hours after plating, the medium was replaced with fresh medium supplemented by the PPARγ ligands troglitazone (10 μM) (Biomol, Plymouth Meeting, PA) or GW7845 (1 μM, a gift from GlaxoWellcome), or by dimethylsulfoxide (DMSO) vehicle control. Fresh media and ligands were added after 24 h in culture. After 48 h the cells were collected for RNA. Total RNA was isolated using Tri-reagent (MRC, Cincinnati, OH, USA) and purified using RNeasy (Qiagen, Valencia, CA, USA). RNA samples (30 μg) from three placentas were mixed, and the mixture was used for double stranded cDNA synthesis using Superscript Choice system (InVitrogen Life Technologies, Carlsbad, CA, USA) and a T7T21 oligonucleotide primer (GenSet, La Jolla, CA, USA). Biotin-labeled RNA was synthesized by in vitro transcription using Enzo Bioarray RNA labeling kit (Enzo Diagnostics, Farmingdale, NY, USA). The RNA was fragmented and divided into five identical aliquots (15 μg of cRNA per aliquot). Each of the five cRNA aliquots was added to hybridization cocktail and sequentially hybridized to identical lot number U95A, U95B, U95C, U95D, and U95E arrays. All arrays were hybridized, washed, stained, and scanned using standard Affymetrix protocols. Together, we used identical cRNA samples to probe the expression of approximately 60 000 genes for each of the three conditions, for a total of 180 000 probe sets, each sampled in five replicates using a total of 75 chips.
Standard Affymetrix protocols using control oligonucleotide B2 were used for proper scanning and grid alignment, and samples were pre-tested using a control cRNA mix from Escherichia coli bioB, bioC, bioD, and P1 cre recombinase, used for monitoring of hybridization, washing, and staining conditions as well as reference samples for normalizing between experiments. Immediately after hybridization the chips were placed in the Affymetrix GeneChip Fluidics Station 400 and sequentially processed for low stringency wash, followed by high stringency wash, streptavidin / phycoerythrin stain, repeat low stringency wash, anti-streptavidin antibody stain, a second streptavidin / phycoerythrin stain and a final low stringency wash. After washing and staining each chip was placed in the Affymetrix Gene-Chip array scanner for image capture and conversion to numerical output using the Microarray Analysis Suite version 5.0. Comparison between chips using the Affymetrix protocol was performed using 'baseline chip' intensity values, normalized to average signal intensity.
This research was supported by NIH R01 ES11597-01 and the Siteman Cancer Center GeneChip Core Facility, Washington University School of Medicine, St Louis, MO, USA. We thank Elena Sadovsky and Lori Rideout for technical assistance. We thank Tim Willson from GlaxoWellcome for generously providing GW7845. We also thank Jeff Milbrandt and Mark Watson (both from Washington University) for discussions.
- Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA. 2001, 98: 5116-5121. 10.1073/pnas.091062498.PubMed CentralView ArticlePubMedGoogle Scholar
- Yang YH, Dudoit S, Luu P, Lin DM, Peng V, Ngai J, Speed TP: Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucleic Acids Res. 2002, 30: 1-10.View ArticleGoogle Scholar
- Tseng GC, Oh MK, Rohlin L, Liao JC, Wong WH: Issues in cDNA microarray analysis: quality filtering, channel normalization, models of variations and assessment of gene effects. Nucleic Acids Res. 2001, 29: 2549-2557. 10.1093/nar/29.12.2549.PubMed CentralView ArticlePubMedGoogle Scholar
- Li C, Wong WH: Model-based analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Natl Acad Sci USA. 2001, 98: 31-36. 10.1073/pnas.011404098.PubMed CentralView ArticlePubMedGoogle Scholar
- Wolfinger RD, Gibson G, Wolfinger ED, Bennett L, Hamadeh H, Bushel P, Afshari C, Paules RS: Assessing gene significance from cDNA microarray expression data via mixed models. J Comput Biol. 2001, 8: 625-637. 10.1089/106652701753307520.View ArticlePubMedGoogle Scholar
- Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP: Summaries of Affymetrix GeneChip probe level data. Nucleic Acids Res. 2003, 31: e15-10.1093/nar/gng015.PubMed CentralView ArticlePubMedGoogle Scholar
- Duggan DJ, Bittner M, Chen Y, Meltzer P, Trent JM: Expression profiling using cDNA microarrays. Nat Genet. 1999, 21: 10-14. 10.1038/4434.View ArticlePubMedGoogle Scholar
- Kerr MK, Martin M, Churchill GA: Analysis of variance for gene expression microarray data. J Comput Biol. 2000, 7: 819-837. 10.1089/10665270050514954.View ArticlePubMedGoogle Scholar
- Satagopan JM, Panageas KS: A statistical perspective on gene expression data analysis. Stat Med. 2003, 22: 481-499. 10.1002/sim.1350.View ArticlePubMedGoogle Scholar
- Churchill GA: Fundamentals of experimental design for cDNA microarrays. Nat Genet. 2002, 32 Suppl: 490-495. 10.1038/ng1031.View ArticlePubMedGoogle Scholar
- Chen Y, Dougherty ER, Bittner ML: Ratio-based decisions and the quantitative analysis of cDNAmicroarray images. J Biomed Optics. 1997, 2: 364-374. 10.1117/1.429838.View ArticleGoogle Scholar
- Mariani TJ, Budhraja V, Mecham BH, Gu CC, Watson MA, Sadovsky Y: A variable fold change threshold determines significance for expression microarrays. FASEB J. 2003, 17: 321-323.PubMedGoogle Scholar
- Tsien CL, Libermann TA, Gu X, Kohane IS: On reporting fold differences. Pac Symp Biocomput. 2001, 496-507.Google Scholar
- Durbin BP, Hardin JS, Hawkins DM, Rocke DM: A variance-stabilizing transformation for gene expression microarray data. Bioinformatics. 2002, 18: S105-110.View ArticlePubMedGoogle Scholar
- Baldi P, Long AD: A Bayesian framework for the analysis of microarray expression data: regularized t -test and statistical inferences of gene changes. Bioinformatics. 2001, 17: 509-519. 10.1093/bioinformatics/17.6.509.View ArticlePubMedGoogle Scholar
- Efron B, Tibshirani R, Storey JD, Tusher VG: Empirical Bayes analysis of a microarray experiment. J Am Statistical Assoc. 2001, 96: 1151-1160. 10.1198/016214501753382129.View ArticleGoogle Scholar
- Cui X, Churchill GA: Statistical tests for differential expression in cDNA microarray experiments. Genome Biol. 2003, 4: 210-10.1186/gb-2003-4-4-210.PubMed CentralView ArticlePubMedGoogle Scholar
- Cleveland W, Devlin S: Locally-weighted regression: an approach to regression analysis by local fitting. J Am Statistical Assoc. 1988, 83: 596-610.View ArticleGoogle Scholar
- Colantuoni C, Henry G, Zeger S, Pevsner J: Local mean normalization of microarray element signal intensities across an array surface: quality control and correction of spatially systematic artifacts. Biotechniques. 2002, 32: 1316-1320.PubMedGoogle Scholar
- Edwards D: Non-linear normalization and background correction in one-channel cDNA microarray studies. Bioinformatics. 2003, 19: 825-833. 10.1093/bioinformatics/btg083.View ArticlePubMedGoogle Scholar
- Rocke DM, Durbin BP: A model for measurement error for gene expression arrays. J Comput Biol. 2001, 8: 557-569. 10.1089/106652701753307485.View ArticlePubMedGoogle Scholar
- Lee ML, Kuo FC, Whitmore GA, Sklar J: Importance of replication in microarray gene expression studies: statistical methods and evidence from repetitive cDNA hybridizations. Proc Natl Acad Sci USA. 2000, 97: 9834-9839. 10.1073/pnas.97.18.9834.PubMed CentralView ArticlePubMedGoogle Scholar
- Newton MA, Kendziorski CM, Richmond CS, Blattner FR, Tsui KW: On differential variability of expression ratios: improving statistical inference about gene expression changes from microarray data. J Comput Biol. 2001, 8: 37-52. 10.1089/106652701300099074.View ArticlePubMedGoogle Scholar
- Ideker T, Thorsson V, Siegel AF, Hood LE: Testing for differentially-expressed genes by maximum-likelihood analysis of microarray data. J Comput Biol. 2000, 7: 805-817. 10.1089/10665270050514945.View ArticlePubMedGoogle Scholar
- Westfall PH, Young SS: p-value adjustment for multiple tests in multivariate binomial models. J Am Statistical Assoc. 1989, 84: 780-786.Google Scholar
- Reiner A, Yekutieli D, Benjamini Y: Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003, 19: 368-375. 10.1093/bioinformatics/btf877.View ArticlePubMedGoogle Scholar
- Kliman HJ, Nestler JE, Sermasi E, Sanger JM, Strauss JM: Purification, characterization and in vitro differentiation of cytotrophoblasts from human term placentae. Endocrinology. 1986, 118: 1567-1582.View ArticlePubMedGoogle Scholar
- Schaiff WT, Carlson MG, Smith SD, Levy R, Nelson DM, Sadovsky Y: Peroxisome proliferator-activated receptor-g modulates differentiation of human trophoblast in a ligand-specific manner. J Clin Endocrinol Metab. 2000, 85: 3874-3881. 10.1210/jc.85.10.3874.PubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.