Statistical modeling for selecting housekeeper genes
© Szabo et al.; licensee BioMed Central Ltd. 2004
Received: 16 January 2004
Accepted: 11 June 2004
Published: 29 July 2004
There is a need for statistical methods to identify genes that have minimal variation in expression across a variety of experimental conditions. These 'housekeeper' genes are widely employed as controls for quantification of test genes using gel analysis and real-time RT-PCR. Using real-time quantitative RT-PCR, we analyzed 80 primary breast tumors for variation in expression of six putative housekeeper genes (MRPL19 (mitochondrial ribosomal protein L19), PSMC4 (proteasome (prosome, macropain) 26S subunit, ATPase, 4), SF3A1 (splicing factor 3a, subunit 1, 120 kDa), PUM1 (pumilio homolog 1 (Drosophila)), ACTB (actin, beta) and GAPD (glyceraldehyde-3-phosphate dehydrogenase)). We present appropriate models for selecting the best housekeepers to normalize quantitative data within a given tissue type (for example, breast cancer) and across different types of tissue samples.
Genes that exhibit minimal variation in messenger RNA (mRNA) quantity across a variety of cell types and biological conditions provide valuable controls for relative quantification. Normalizing quantitative data with housekeeper genes has many applications, from identifying genes regulated during embryogenesis to developing new cancer diagnostics. Although finding biological significance in gene-expression data can rely heavily on the performance of the housekeeper genes, there is a paucity of information on testing these genes for their suitability for this role.
The copy number of a housekeeper gene should be proportional to the amount of poly(A) RNA present in the sample and this proportion should be maintained across a variety of experimental conditions. As nucleic acids show high absorbance at 260 nm (A260), spectrophotometers provide approximate amounts of total DNA/RNA present in a sample. Using absorbance methods alone, however, gives no information about the type of nucleic acid (for example DNA versus RNA) or contributions from different nucleic acid fractions (for example rRNA versus mRNA). It is assumed that mRNA comprises approximately 1-3% of the total RNA. However, this contribution may change depending on the extraction method used. For instance, column extraction methods provide better exclusion of ribosomal RNA than solvent extraction methods . By combining capillary electrophoresis with absorbance, it is possible to accurately quantify these different fractions .
Traditionally, housekeepers have been used in Northern blot analysis to represent the amount of mRNA in the sample and to control for sample loading, blot transfer and probe hybridization. Highly expressed genes serving fundamental roles in the cell are commonly used for this purpose but may not be optimal under certain experimental conditions [3–5]. For example, the sensitivity and accuracy of northern blot analysis with densitometry may be decreased using a highly expressed housekeeper gene that can saturate the autoradiographic signal . To resolve this problem and compensate for limitations in dynamic range, control genes may be chosen to have a level of gene expression similar to the gene(s) of interest (that is, the test genes).
Microarrays are more practical for genome-wide expression analysis than northern blots . With cDNA microarrays, a common reference sample is usually used to compare the expression of each gene across many experimental sample(s) [8, 9]. Because each gene in the experimental sample is directly compared to the same gene in the common reference, housekeeper genes are not necessary for normalization. Microarrays are commonly applied to finding genes with differential expression across experimental conditions, however the data may also be used to identify stably expressed genes that can serve as important controls for northern blot analysis, ribonuclease protection assays and quantitative reverse transcription PCR (RT-PCR). In turn, these other quantitative methods are often used to verify differentially expressed genes identified by microarray [10–12].
Housekeeper genes are often adopted from the literature and used across a variety of experimental conditions, some of which may induce differences in their expression. If unrecognized, unexpected changes in housekeeper expression could result in erroneous conclusions about real biological effects such as responses to drugs. In addition, this type of change would be difficult to detect because most experiments only include a single housekeeper gene. It is difficult to determine whether a given gene has the constitutive property of a housekeeper when the true amount of mRNA in a sample is unknown. As a way round this dilemma, Vandesompele et al. postulate that gene pairs that have stable expression patterns relative to each other are proper control genes . An alternative method for quantitative analysis of RT-PCR data that does not require housekeeper genes for normalization is to use global pattern recognition (GPR). For instance, Akilesh et al. used a GPR algorithm to search for eligible normalizing genes within an assay plate and then used those genes as controls to identify differentially expressed genes .
Although relative quantification using housekeeper genes is a practical method of estimating the expression level of a test gene, the transcript amount in the sample is a summation and the method does not consider transcript differences on a cell-to-cell basis. Fluorescence in situ hybridization (FISH) is clinically used to determine absolute DNA copy number (for example, HER2 amplification) in a cell, but these methods still average the copy number after counting many cells and the technique is expensive and laborious . In situ methods for detecting RNA transcripts have been developed but the assays are semiquantitative and subjective .
In the work presented here, we applied several models to selecting the best housekeeper genes for breast cancer and give algorithms that can be generalized to find housekeeper genes that are appropriate for normalizing quantitative data within and between tissue types.
Results and discussion
One tissue type
To select the best housekeepers for normalizing data across a single tissue type, we tested three variations of a model (Model 1, a-c) with real-time quantitative RT-PCR data generated from primary breast samples (see Materials and methods for details).
We model the expression y ij of gene j in sample I by
where μ is the overall mean (log-) expression, T i is the difference of the ith tissue sample from the overall average and G j is the difference of the jth gene from the overall average. The key feature of this model that makes it different from a traditional ANOVA model is that it allows heteroscedastic errors to account for different variability in the genes . The variability around the gene-specific mean log-expression μ + T i + G j is quantified by the error standard deviation σ j . The Bayesian information criterion (BIC) was used to avoid overfitting the data . Model 1a had the best BIC value and was selected from a range of competing models that included a method with equal error variances (Model 1b in Materials and methods) and a more complex method with correlated errors (Model 1c in Materials and methods).
Standard deviation estimates of log expression using Model 1a for selecting the single best housekeeper gene for breast cancer
Estimated standard deviation
95% confidence interval
Standard deviation estimates of log expression using Model 1a for selecting the best housekeeper gene(s) for breast cancer
PSMC4, MRPL19, PUM1
PSMC4, MRPL19, PUM1, SP3A1
PSMC4, MRPL19, SF3A1, PUM1, ACTB
PSMC4, MRPL19, SF3A1, PUM1, GAPD, ACTB
These findings illustrate the importance of performing an unbiased and genome-wide search for housekeepers rather than relying on traditional housekeeper genes. We used microarray data to select genes with low variability in expression across breast tumors and cell lines. Because the quantitative differences between the microarray and RT-PCR platforms are relative, genes with low variability in expression across tumors by microarray should also show low variability in expression by RT-PCR. Although the quantitative data from microarray tends to have an overall smaller dynamic range compared to RT-PCR, this is primarily due to loss of information from genes expressed at low levels. Our microarray dataset was filtered to remove genes with signals near background noise.
The result is very similar using Vandesompele et al.'s M value method, with only the positions of PUM1 and PSMC4 changing in stability rank. It should be noted that the M-value method does not order the two best genes (MRPL19 and PSMC4). Their best gene-set selection approach would suggest using the (log-scale) average of these two best genes as a control. Such a concordance is not surprising given the close relationship between the M value and our model using the variability of the average of several genes (see Materials and methods for details). A benefit of our approach is the ability to compare the variability of individual genes to that of an average of several genes.
Multiple tissue types
Gene(s) with minimal variation in expression across different cell types serve as good 'universal' housekeepers. A universal control may be a single gene or combination of genes. While the former should display both low variability within a given tissue type and consistent basal levels of expression across tissue types, the latter may comprise a gene set with individually different, but complementary, basal expression levels across tissue types.
To compare the performance of housekeepers within and between different tissues, we made a Model 2 (see Materials and methods for further details) that models the expression of gene j in the ith sample of tissue-type k by
Components of the standard deviation estimates of the log-expression of the data of Vandesompele et al. 
Standard deviation of genes (σ j )
Tissue-type specific multipliers (ς k )
Cultured normal fibroblasts
Neuroblastoma cell lines
Minimax MSE optimal gene sets for each set size
Maximum number of members
HPRT1, RPL13A, UBC
HPRT1, RPL13A, UBC
HPRT1, RPL13A, UBC
ACTB, HPRT1, SDHA, TBP, UBC, YWHAZ
ACTB, HPRT1, RPL13A, SDHA, TBP, UBC, YWHAZ
ACTB, HPRT1, RPL13A, SDHA, TBP, UBC, YWHAZ
ACTB, HPRT1, RPL13A, SDHA, TBP, UBC, YWHAZ
ACTB, B2M, GAPD, HMBS, HPRT1, RPL13A, SDHA, TBP, UBC, YWHAZ
In summary, we have modeled the performance of putative housekeepers to test their goodness-of-fit in serving as normalization controls for relative insert quantification. A major advantage of a model approach is that the terms are placed within a solid statistical framework and are not ad hoc, which allows the algorithm to be generalized to a variety of different experimental conditions. The genes and algorithms that we have selected for normalization should have broad utility for diagnostics and research.
Materials and methods
Pre-selection of assayed genes from microarray experiments
Four candidate housekeepers (PSMC4, MRPL19, PUM1 and SF3A1) were selected from a microarray dataset containing 40 different breast tumors, three normal breast samples and 19 cell lines representing 17 different cell lines of diverse nature including lymphocytes, fibroblasts and epithelial cells . All experiments were done using a common reference strategy in which all experimental samples are compared to the same reference comprised of a pool of RNAs isolated from 11 diverse human cell lines .
To select housekeepers, we first 'filtered' the microarray data to select genes with Cy3 and Cy5 signal intensities greater than 500 units across at least 75% of the experiments. This requirement ensures that the gene is well expressed not only in the experimental samples, but also in the common reference sample. Next, we used the SAS/STAT Analysis Package Version 8 (SAS Institute Inc., Cary, NC) to identify a set of genes that showed a small range of expression across sample types and the least variance of the array-mean normalized log-ratios. For real-time RT-PCR, we selected four of the top six genes - PUM1, PSMC4, MRPL19 and SF3A. The two other low-variability genes identified in the data were IER3 (immediate early response 3) and SRY ((sex determining region Y)-box 2). We did not select these genes because of their potential for being differentially regulated under other conditions. However, we did include GAPD and ACTB, which are commonly used reference genes , in the set of candidate genes for comparison to the microarray selection.
Samples and cDNA preparation
Breast samples were acquired under informed consent and received at the Huntsman Cancer Institute (Salt Lake City, UT) for gene expression analysis (University of Utah, IRB #8533). All specimens were expediently processed in pathology upon arrival from surgery. Samples were grossly dissected, procured by flash freezing in liquid nitrogen, and stored at -80°C until RNA extraction. Approximately 50-100 mg cancer tissue was homogenized from each sample, and total RNA was prepared using the RNeasy midi kit (Qiagen). The integrity of RNA was determined using the RNA 6000 Nano LabChip kit (Agilent Technologies) and an Agilent 2100 Bioanalyzer. Two microliters of total RNA (50 ng/μl) were heated to 70°C and 1 μl was loaded on the column. Degradation was evaluated using the signal of the 18S and 28S ribosomal peaks .
First-strand cDNA was synthesized from 1 μg total RNA using oligo(dT) primers and Superscript III reverse transcriptase following manufacturer's instructions (Superscript III First-Strand Synthesis System, Invitrogen Life Technologies). Briefly, the reaction was held at 48°C for 50 min, followed by a 15 min step at 70°C. The cDNA was washed on QIAquick PCR purification column (Qiagen) and eluted in 2 × 50 μl of elution buffer. The cDNA was then diluted in TE' (10 mM Tris, 0.1 mM EDTA, pH 8.0), aliquoted and stored at -80°C for further use.
Real-time quantitative PCR
Primers for housekeeper genes
PSMC4 (UniGene reference Hs.211594 - Gene ID: 5704)
MRPL19 (UniGene reference Hs.44024 - Gene ID: 9801)
PUM1 (UniGene reference Hs.153834 - Gene ID:9698)
SF3A1 (UniGene reference Hs.406277 -Gene ID:10291)
ACTB (UniGene reference Hs.426930 - Gene ID:60)
†GAPD (UniGene reference Hs.169476 - Gene ID:2597)
PCR was done using the following protocol: initial denaturation 95°C for 1 min 30 sec, then 50 cycles at 94°C for 1 sec for denaturation, 60°C for 5 sec (20°C/sec transition) for annealing, 72°C for 8 sec (2°C/sec transition) for extension. Fluorescence emission of SYBR Green I (channel 1, 530 nm) was acquired each cycle after the extension step. A melting step was performed after PCR to determine product purity. For melting curve analysis, the reactions were rapidly (20°C/sec) cooled from 95°C to 60°C and then slowly heated (0.1°C/sec) back to 95°C while continuously monitoring fluorescence.
Copy number was determined using the crossing point (Cp) value, which is automatically calculated using the LightCycler 3.5 software (Roche Molecular Biochemicals). The Cp value is reported as a fractional cycle number that is determined from the second derivative maximum (point of maximum acceleration) on the PCR amplification curve (fluorescence versus cycle number) . A relative starting copy number was determined for each housekeeper using a calibration curve done with the same batch of master mix. Efficiency (E) of PCR was calculated from a plot of Cp versus log ng cDNA .
E = 10-1/slope
Modeling expression data
As the effects of interest are fold changes, we modeled the log-transformed expression Model 1a.
log y ij = μ + T i + G j + ε ij ,
where μ denotes the overall mean (log) expression, T i is the difference of the ith tissue sample from the overall average and G j is the difference of the jth gene from the overall average. The key feature of this model that makes it different from a traditional ANOVA model is that it allows heteroscedastic errors: the variability of the genes is different.
We fitted the model using the gls routine of the nlme library for R, however other commonly available software such as PROC MIXED from SAS could have been used.
Based on the model, the variability of the logarithm of the geometric mean
of a gene-set S was estimated as
Vandesompele et al.'s M-value is the average of relative standard deviations of the log-expression levels. Under Model 1, the M-value of the gene is closely related to its variance (under Models 2 and 3 below, the similar relationships can be derived):
We tested the assumption of unequal variances by fitting Model 1b that forces all the genes to have the same variability (this is the classical ANOVA model).
log y ij = μ + T i + G j + ε ij ,
Model 1c with a correlated error structure can be used to assess the assumption of (conditional) independence of the genes given the sample mean. If warranted, a more complicated correlation structure can be imposed.
log y ij = μ + T i + G j + ε ij ,
For the multiple tissue-type set-up the notation and the model need to be extended. We will denote the expression level of gene j of in the ith sample of type k by yi(k)j, i = 1,...n k , j = 1,...,g and k = 1,...,m. The best-fitting model for the data, which we call Model 2, had the form
log y i(k)j = μ + C k + T i(k) + G j + (CG) kj + ε i(k)j ,
Thus the errors are independent and their variability is decomposed into a gene-specific and tissue-type specific multiplicative components. The last restriction ensures the uniqueness of the solution. Simpler models that we considered used uniform error variance, equal error variance for tissue types, and equal error variance for genes. We also considered more complex models that used exchangeable correlation structure for the errors and unstructured error variance (each gene-tissue-type combination has a variance parameter). The BIC was used as a basis for model selection.
Additional data files
Additional data available with this paper online is an Excel file with the relative copy numbers of six genes in the 80 breast cancer samples used in this study (Additional data file 1).
This work has been supported in part by the National Cancer Institute (R33 CA097769-01).
- Miller CL, Yolken RH: Methods to optimize the generation of cDNA from postmortem human brain tissue. Brain Res Brain Res Protoc. 2003, 10: 156-167. 10.1016/S1385-299X(02)00214-3.PubMedView ArticleGoogle Scholar
- Panaro NJ, Yuen PK, Sakazume T, Fortina P, Kricka LJ, Wilding P: Evaluation of DNA fragment sizing and quantification by the Agilent 2100 bioanalyzer. Clin Chem. 2000, 46: 1851-1853.PubMedGoogle Scholar
- Suzuki T, Higgins PJ, Crawford DR: Control selection for RNA quantitation. Biotechniques. 2000, 29: 332-337.PubMedGoogle Scholar
- Bhatia P, Taylor WR, Greenberg AH, Wright JA: Comparison of glyceraldehyde-3-phosphate dehydrogenase and 28S-ribosomal RNA gene expression as RNA loading controls for northern blot analysis of cell lines of varying malignant potential. Anal Biochem. 1994, 216: 223-226. 10.1006/abio.1994.1028.PubMedView ArticleGoogle Scholar
- Spanakis E: Problems related to the interpretation of autoradiographic data on gene expression using common constitutive transcripts as controls. Nucleic Acids Res. 1993, 21: 3809-3819.PubMedPubMed CentralView ArticleGoogle Scholar
- Eggert A, Brodeur GM, Ikegaki N: Relative quantitative RT-PCR protocol for TrkB expression in neuroblastoma using GAPD as an internal control. Biotechniques. 2000, 28: 681-691.PubMedGoogle Scholar
- Schena M, Shalon D, Davis RW, Brown PO: Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science. 1995, 270: 467-470.PubMedView ArticleGoogle Scholar
- Perou CM, Sorlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, Pollack JR, Ross DT, Johnsen H, Akslen LA, et al: Molecular portraits of human breast tumours. Nature. 2000, 406: 747-752. 10.1038/35021093.PubMedView ArticleGoogle Scholar
- van de Vijver MJ, He YD, van't Veer LJ, Dai H, Hart AA, Voskuil DW, Schreiber GJ, Peterse JL, Roberts C, Marton MJ, et al: A gene-expression signature as a predictor of survival in breast cancer. N Engl J Med. 2002, 347: 1999-2009. 10.1056/NEJMoa021967.PubMedView ArticleGoogle Scholar
- Dhanasekaran SM, Barrette TR, Ghosh D, Shah R, Varambally S, Kurachi K, Pienta KJ, Rubin MA, Chinnaiyan AM: Delineation of prognostic biomarkers in prostate cancer. Nature. 2001, 412: 822-826. 10.1038/35090585.PubMedView ArticleGoogle Scholar
- Welsh JB, Zarrinkar PP, Sapinoso LM, Kern SG, Behling CA, Monk BJ, Lockhart DJ, Burger RA, Hampton GM: Analysis of gene expression profiles in normal and neoplastic ovarian tissue samples identifies candidate molecular markers of epithelial ovarian cancer. Proc Natl Acad Sci USA. 2001, 98: 1176-1181. 10.1073/pnas.98.3.1176.PubMedPubMed CentralView ArticleGoogle Scholar
- Mischel PS, Nelson SF, Cloughesy TF: Molecular analysis of glioblastoma: pathway profiling and its implications for patient therapy. Cancer Biol Ther. 2003, 2: 242-247.PubMedView ArticleGoogle Scholar
- Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3: research0034.1-0034.11. 10.1186/gb-2002-3-7-research0034.View ArticleGoogle Scholar
- Akilesh S, Shaffer DJ, Roopenian D: Customized molecular phenotyping by quantitative gene expression and pattern recognition analysis. Genome Res. 2003, 13: 1719-1727. 10.1101/gr.533003.PubMedPubMed CentralView ArticleGoogle Scholar
- Tubbs RR, Pettay JD, Roche PC, Stoler MH, Jenkins RB, Grogan TM: Discrepancies in clinical laboratory testing of eligibility for trastuzumab therapy: apparent immunohistochemical false-positives do not get the message. J Clin Oncol. 2001, 19: 2714-2721.PubMedGoogle Scholar
- Kristt D, Turner I, Koren R, Ramadan E, Gal R: Overexpression of cyclin D1 mRNA in colorectal carcinomas and relationship to clinicopathological features: an in situ hybridization analysis. Pathol Oncol Res. 2000, 6: 65-70.PubMedView ArticleGoogle Scholar
- Pinheiro JCBD: Mixed-effects Models in S and S-PLUS. 2000, New York: SpringerView ArticleGoogle Scholar
- Schwarz G: Estimating the dimension of a model. Annls Stat. 1978, 6: 461-464.View ArticleGoogle Scholar
- Perou CM, Brown PO, Botstein D: Tumor classification using gene expression patterns from DNA microarrays. New Technologies for Life Sciences: A Trends Guide. 2000, 67-76.Google Scholar
- Roux S, Pichaud F, Quinn J, Lalande A, Morieux C, Jullienne A, de Vernejoul MC: Effects of prostaglandins on human hematopoietic osteoclast precursors. Endocrinology. 1997, 138: 1476-1482. 10.1210/en.138.4.1476.PubMedView ArticleGoogle Scholar
- Frank SG, Bernard PS: Profiling breast cancer using real-time quantitative PCR. In Rapid Cycle Real-Time PCR: Methods and Applications. Edited by: Wittwer CT, Meuer S, Nakagawara K. 2003, Heidelberg: Springer, 95-106.Google Scholar
- Rasmussen RP: Quantification on the LightCycler. In Rapid Cycle Real-Time PCR: Methods and Applications. Edited by: Wittwer CT, Meuer S, Nakagawara K. 2003, Heidelberg: Springer, 21-34.Google Scholar
- SantaLucia J: A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. Proc Natl Acad Sci USA. 1998, 95: 1460-1465. 10.1073/pnas.95.4.1460.PubMedPubMed CentralView ArticleGoogle 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.