|Home | About | Journals | Submit | Contact Us | Français|
An important consideration when analyzing both microarray and quantitative PCR expression data is the selection of appropriate genes as endogenous controls or reference genes. This step is especially critical when identifying genes differentially expressed between datasets. Moreover, reference genes suitable in one context (e.g. lung cancer) may not be suitable in another (e.g. breast cancer). Currently, the main approach to identify reference genes involves the mining of expression microarray data for highly expressed and relatively constant transcripts across a sample set. A caveat here is the requirement for transcript normalization prior to analysis, and measurements obtained are relative, not absolute. Alternatively, as sequencing-based technologies provide digital quantitative output, absolute quantification ensues, and reference gene identification becomes more accurate.
Serial analysis of gene expression (SAGE) profiles of non-malignant and malignant lung samples were compared using a permutation test to identify the most stably expressed genes across all samples. Subsequently, the specificity of the reference genes was evaluated across multiple tissue types, their constancy of expression was assessed using quantitative RT-PCR (qPCR), and their impact on differential expression analysis of microarray data was evaluated.
We show that (i) conventional references genes such as ACTB and GAPDH are highly variable between cancerous and non-cancerous samples, (ii) reference genes identified for lung cancer do not perform well for other cancer types (breast and brain), (iii) reference genes identified through SAGE show low variability using qPCR in a different cohort of samples, and (iv) normalization of a lung cancer gene expression microarray dataset with or without our reference genes, yields different results for differential gene expression and subsequent analyses. Specifically, key established pathways in lung cancer exhibit higher statistical significance using a dataset normalized with our reference genes relative to normalization without using our reference genes.
Our analyses found NDUFA1, RPL19, RAB5C, and RPS18 to occupy the top ranking positions among 15 suitable reference genes optimal for normalization of lung tissue expression data. Significantly, the approach used in this study can be applied to data generated using new generation sequencing platforms for the identification of reference genes optimal within diverse contexts.
Gene expression profiling, including quantitative RT-PCR (qPCR) and microarray experimentation, is invaluable for the molecular analysis of biological systems. The interpretation of results from such experiments (i.e., the determination of differential expression for a particular gene among datasets) is strongly influenced by the selection of reference genes for normalization across datasets . Specifically, gene expression is normalized within a given dataset by calculating the transcript abundance of the gene of interest relative to a gene that is constantly expressed across independent datasets (termed a "housekeeping" or a "reference" gene), and differential expression between two datasets or samples is determined by calculating the ratio of the normalized expression levels for the gene of interest between the two datasets. Typically, housekeeping genes satisfy the following criteria: they are highly expressed in the cell, the variability in expression between samples is minimal, and the genes' expression is not influenced by the experimental conditions tested . Hence, problems arise when housekeeping genes are selected that do not meet these criteria, as fluctuations in these genes may erroneously influence the data interpretation.
Historically, beta actin (ACTB), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and 18 S rRNA have been routinely used as reference genes for qPCR and microarray data normalization. However, a number of studies have shown that expression of these genes varies considerably depending on the specific tissue type and disease state of the tissue [3-16]. Attempts to achieve more reliable normalization include the spiking of synthetic poly-A RNAs for the analysis of cDNA arrays and northern blots, and the combined use of an oligo-(dT)n primer with an 18 S specific primer for qPCR analysis [17,18]. In addition, re-mining of large microarray datasets for the identification of novel, highly stable genes, as well as use of a combination of reference genes instead of a single gene for normalization, are some of the other approaches taken to address this problem [11,13,19].
Recently, efforts have been made to identify more suitable reference genes for microarray and qPCR studies of lung cancer. Specifically, candidate reference genes have been identified from the mining of microarray gene expression data to identify the least variable genes, followed by validation of expression using qPCR [11,20,21]. However, as microarray data do not provide absolute abundance values for transcripts, selection of reference genes from this type of data is inherently problematic. To circumvent this handicap in the utilization of microarray data, we turn to the use of large-scale expression profiling permitted by serial analysis of gene expression (SAGE) experimentation for the identification of novel reference genes optimal for the study of lung cancer. This approach, which we have termed normalization of expression by permutation of SAGE (NEPS), takes advantage of the fact that SAGE is a transcriptome profiling technique that identifies the absolute abundance levels of transcripts by direct enumeration of sequence tag counts, thus allowing the direct comparison of expression levels across multiple profiles without the need for reference or housekeeping genes .
NEPS adopts a permutation test approach designed for analyzing relatively small sample sizes, such as those typically encountered with SAGE. Unlike the conventional T-test, the permutation test is non-parametric . The null hypothesis states that the mean gene expression levels in two groups of SAGE libraries being compared (in this case normal and cancer), are the same. For this analysis, samples from both the normal and the cancer groups are pooled, followed by random sampling to create a simulated Group 1 and a simulated Group 2. For each gene, the difference in expression between these two simulated groups was measured. This exercise was repeated 10,000 times, thus generating a simulated mean μ and a simulated standard deviation σ. The permutation score (PS) of a given gene is defined by , where O is the true difference between the average expression levels in the two groups. Hence, for a given gene, the closer the permutation score is to zero, the more it satisfies the constancy requirement.
To demonstrate the utility of NEPS for selecting genes that satisfy the constancy requirement, we analyzed 24 bronchial epithelial lung SAGE libraries, 2 lung parenchyma libraries, and 11 lung squamous cell carcinoma libraries. From this analysis, NEPS selected 15 genes, which we hereafter refer to as the lung-NEPS reference genes (Table (Table1).1). We further demonstrate that (1) while these genes perform well as reference genes for lung, they are not satisfactory for normalization of expression data from other tissues, suggesting that reference genes are tissue-specific, and (2) in lung cancer datasets, differential gene expression determination and subsequent pathway analyses are improved after normalization using the lung-NEPS reference genes.
26 normal and 11 lung cancer SAGE libraries were constructed and used in the analysis . The construction of the 26 normal libraries has been previously described [25,26]. 24 of these libraries were generated from exfoliated bronchial cells obtained from bronchial brushes, and two libraries from normal lung parenchyma (Additional file 1). Conversely, the 11 cancer libraries were generated from biopsied specimens with six libraries representing lung squamous cell carcinoma and five libraries representing carcinoma in situ. This data can be found at the GEO database with the following series accession numbers: GSE3707, GSE5473, and GSE7898. All samples were acquired under approval by the University of British Columbia - British Columbia Cancer Agency Research Ethics Board (UBC-BCCA-REB) and all subjects provided written consent.
Publicly available SAGE data were also used in this analysis, representing both brain and breast cancer. Specifically, six normal and 12 breast cancer libraries (Additional file 2) and 7 normal and 19 brain cancer libraries were used (Additional file 3). The libraries were obtained from the cancer genome anatomy project (CGAP) database http://cgap.nci.nih.gov[27,28].
Given that SAGE libraries are expensive to generate, the number of libraries in a given study is typically small (i.e., in 10's, rather than in 100's). The permutation test is a non-parametric test, which does not assume any underlying distribution. The number of samples required for the test to achieve sufficient statistical power is relatively low compared to other statistical tests (e.g., t-test and χ2-square test). Furthermore, each additional sample increases the power of the test exponentially. The permutation test is a test of the means between two different distributions. Without loss of generality, let us assume that one distribution is for the gene expression level of a particular gene in normal tissues (i.e. subscript n), and that the other distribution is for cancerous tissues (i.e. subscript c). Genes are selected using the following hypotheses:
If there is little difference between the two means, it would make no difference if we mix the cancerous samples with the normal samples. But, if the null hypothesis is rejected, it indicates that the gene expression levels of normal and cancer samples are sufficiently different (the alternative hypothesis). In the following, we show our specific implementation of the test. Let n and c be the number of normal tissue samples and the number of cancerous tissue samples respectively.
A. For each gene, select all the gene-specific normalized tag counts from the normal libraries and all the gene-specific normalized tag counts from the cancer libraries.
B. Randomly select n counts to create a simulated normal set, and calculate the simulated normal mean μsn.
C. Similarly, select the remaining c counts form the simulated cancerous set. Calculate the simulated cancer mean μsc.
D. Consider the random variable v = μsc - μsn, called the simulated difference.
E. Repeat the steps A to D above m times. Let μ and σ denote the mean and the standard deviation of v.
F. Now separate the libraries back into their true identity: normal or cancerous. Calculate the true observed difference O = μrc - μrn, where μrc denotes the true mean count of the cancerous libraries, and μrn denotes the true mean of the normal libraries.
G. Calculate the Permutation Score PS where .
H. Repeat all the above steps for each gene. Sort the permutation score in descending order.
The permutation score is one way to measure how likely the actual observed difference occurs by chance. It is based on standardization, i.e., subtracting the mean and then divided by the standard deviation. The more the true observed difference is from the average (expressed as multiples of the standard deviation), the less likely that the true observed difference is a coincidence. That is to say, the larger the permutation score, the more significant is the observed difference between cancerous and normal samples.
On the other hand, for the sake of evaluating the constancy requirement, the ideal reference gene would have a permutation score equal to 0. This means that there is no difference in the distributions of expression levels between cancerous and normal samples. For the results reported here, we used m = 10,000 permutations.
Raw tag counts for each SAGE library were normalized to tags per million (TPM) to facilitate adequate comparison among libraries. Tag-to-gene mapping was performed using the February 5th, 2007 version of SAGEGenie . In cases where multiple SAGE tags mapped to the same gene, the tags were collapsed to capture all potential transcript variants, and a cumulative tag count was utilized for analysis.
The permutation test outlined above was used to identify genes which were statistically similar when comparing the libraries from normal tissue (bronchial epithelium and lung parenchyma) and cancerous tissue of the lung. Three main criteria were used for reference gene selection: permutation score (described above) ≤ 0.15; at least two SAGE tags observed in each library; and an overall average count of ≥ 25 across all samples. For the analysis in brain and breast tissue, the first two criteria were maintained, but due to the lower sequencing depth, an average count of ≥ 10 across all samples was used instead.
One microgram of total RNA from 15 lung tumor and matched non-malignant parenchyma samples were converted to cDNA using the High-Capacity cDNA archive kit (Applied Biosystems Inc., Foster City CA). One hundred nanograms of cDNA were utilized for qPCR using the TaqMan Gene Expression Assay (Applied Biosystems Inc). All fifteen lung NEPS genes and six additional reference genes were assayed. All TaqMan probes were pre-optimized by Applied Biosystems. Primer IDs for all genes are provided in Additional file 4. The 30 samples were assayed in triplicate in parallel along with negative (no cDNA template) controls using the 7500 Fast Real-Time PCR System. Appropriate cDNA dilutions were used such that the exponential phase of the amplification curves were within the 40 PCR cycles recommended by the manufacturer (i.e. ranging from 16-36 cycles for the 20 genes and 1-13 cycles for 18SRNA). Cycle thresholds were determined from amplification curves using 7500 Fast System software.
For the analysis of qPCR data, three different methods were used. Within each method, all genes were ranked from best to worst. Subsequently, for each gene, a cumulative ranking across all three methods was determined by summing its rank from each individual method. Two previously published methods, geNorm  and NormFinder , and the variance of cycle threshold difference (dCt) across all 15 tumor/matched non-malignant sample pairs were the approaches used to determine constancy.
For the Affymetrix data (GSE10072), Raw CEL files were processed through Affymetrix's Microarray Array Suite (MAS) 5.0 algorithm in the "affy" package in Bioconductor [32,33]. Briefly, MAS 5.0 is a three step process which involves a global background signal correction, correction of the probe value for cross-hybridization and spurious signals using mismatch probes which are off by one base, and finally, scale normalization of each experiment to a fixed median intensity to facilitate inter-experimental comparison http://media.affymetrix.com/support/technical/whitepapers/sadd_whitepaper.pdf. Probes were filtered on MAS 5.0 calls, and those having a "P" or "M" call in at least 50% of samples were retained. This resulted in a dataset of 11440 probes. Of the 15 lung NEPS reference genes, 12 were represented on the array platform. Of those 12 genes, probes which had a "P" call in 100% of the samples were used for the calculation of the scaling factor with only one probe/gene allowed. If two probes met these criteria for one gene, the probe with the highest mean expression was chosen. After employing these criteria, eight probes were used (Additional file 5), which represented genes PPP1CB, B2M, RPL4, CAPZB, ATP5J, RAB5C, NDUFA1, and HSPA1A.
For the Agilent microarray data (GSE12428), all lung NEPS genes were represented on this microarray platform. Data was processed as described previously . In the cases where lung NEPS genes were represented with multiple probes, the probe with the maximum average intensity across the dataset was used. A list of the probes used is given in Additional file 6. Since each sample had at least two replicate experiments, the average across replicate experiments was used for each probe.
To determine the scaling factor, for each sample, linear regression analysis was performed comparing the values for the reference gene (x) versus the average values for the reference genes across the sample set (y). The slope of the line based on least-squares fitting was then multiplied to each value in the experiment.
Next, Significance Analysis of Microarrays (SAM) was performed to determine differentially expressed genes between non-malignant and malignant samples for both microarray datasets using the "samr" package in R . Unpaired analysis was performed using the normal samples versus tumor samples and the delta parameter set to 0.4. Probes which had a Q-value% ≤ 5 were considered significant. For the Affymetrix dataset, results were compared between the dataset normalized with MAS 5.0 alone and MAS 5.0 + NEPS scaling and for the Agilent dataset, the comparison was done between median normalization alone and NEPS scaling followed by median normalization.
From our NEPS analysis [with an imposed permutation score (PS) threshold ≤ 0.15, and an average expression of ≥ 25 raw tag counts across all samples], 15 genes were identified as the most consistently expressed across normal and cancerous lung tissue (Table (Table1).1). Here we identified beta-2-microglobulin (B2M), components of the large ribosomal subunit such as ribosomal protein L19 (RPL19) and ribosomal protein L4 (RPL4), components of the small ribosomal subunit such as ribosomal protein S18 (RPS18), and electron transport chain constituents such as NADH dehydrogenase (ubiquinone) 1 alpha subcomplex 1 (NDUFA1), to rate highly in our permutation analysis, thereby suggesting their potential as reliable reference genes. B2 M has previously been utilized as a reference gene [10,15], providing validity to the approach used here. The 18 S and 28 S rRNAs have previously served as reference genes [4,7,11,12,16], and here we show that the ribosomal protein genes can also provide this service.
A previous study reported a meta-analysis of microarray data designed to identify novel reference genes for the study of non-small cell lung cancer (NSCLC) . Using a small panel of tumor/normal specimens, the authors demonstrated that genes commonly used for reference in qPCR experimentation were sub-optimal, and identified novel, more consistently expressed genes to be superior as reference genes. Additionally, studies of asthmatic airways have also shown that traditional reference genes such as ACTB and GAPDH perform poorly in this regard [6,7,12,14,15]. We find similar results using our NEPS-based approach. As described above, genes such as B2 M (permutation score = 0.011) and RPL19 (0.082) were shown to have very low permutation scores denoting stable expression between normal and cancer, whereas ACTB (2.69) and GAPDH (6.48) performed very poorly with significant differential expression. Other known reference genes such as hypoxanthine phosphoribosyltransferase 1 (HPRT1) (0.114) and TATA box binding protein (TBP) (0.468), while exhibiting low permutation scores, were not as highly expressed with average raw tag counts across all samples of 1.54 and 1.65, respectively. In contrast, genes such as peptidylprolyl isomerase A (PPIA, aka cyclophilin A) (6.11), transferrin receptor (TFRC, p90, CD71) (4.70), and phosphoglycerate kinase 1 (PGK1) (5.04) identified in a microarray meta-analysis study (see above) , performed poorly in our study, as revealed by the relatively high permutation scores. Although these particular genes did not perform as well as the reference genes identified from our permutation analysis, other genes identified by Saviozzi et al., such as signal transducer and activator of transcription 1 (STAT1) (0.21), esterase D/formylglutathione hydrolase (ESD) (0.18), Yes-associated protein 1 (YAP1) (0.65) and polymerase (RNA) II (DNA directed) polypeptide A (POLR2A) (0.88) did perform satisfactorily in our study, as evidenced by permutation scores ≤ 1. In addition to these genes, a second set of genes identified using a cross-tissue and cross-platform analysis were also assessed  and similarly, while some genes such as C-terminal binding protein 1 (CTBP1), cullin 1 (CUL1), DIM1 dimethyladenosine transferase 1-like (DIMT1L), tripartite motif-containing 27 (TRIM27) and ubiquilin 1 (UBQLN1) performed reasonably well based on our metric, others such as poly(A) polymerase alpha (PAPOLA) and ADP-ribosylation factor-like 8B (ARL8B) did not (Figure (Figure1,1, Additional file 7).
To further our investigations regarding reference genes optimal for cancer cell biology, we expanded our analysis to include publicly available SAGE libraries representing normal and cancer tissue from both brain and breast. The results of this analysis clearly demonstrate that the reference genes identified in the lung dataset are distinct from those found in either breast (Table (Table2)2) or brain (Table (Table3).3). This data strongly suggests that reference genes should be selected in a tissue specific manner. For example, GAPDH, which performed poorly as a reference gene for lung gene expression analysis (see above), was in fact one of the best reference genes identified from the analysis of the brain dataset (Table (Table3).3). Moreover, not only was there no overlap among the reference gene lists determined for each of the three different tissue types (i.e., lung-NEPS, breast-NEPS, brain-NEPS), but when examining reference genes specific to one tissue type (i.e. lung-NEPS) in the other two tissue types (i.e. breast or brain), the permutation scores for these genes were significantly higher and more variable (Figure (Figure2).2). These results are consistent with other studies demonstrating the need for tissue and context-specific selection of reference genes [3,5,8,10,14].
Using a secondary set of 15 tumor and matched non-malignant samples, qPCR was used to validate consistency of expression for all lung-NEPS genes. Additionally, we performed qPCR for previously identified, commonly used housekeeping genes ACTB, GAPDH, HPRT1, and TBP. In addition, two genes out of 13 identified (CUL1 and TRIM27) as suitable reference genes from a previously published study , were selected here based on high NEPS performance (see above), for qPCR analysis.
Of the NEPS genes analyzed, NDUFA1, RPL19, RAB5C, member RAS oncogene family (RAB5C), and RPS18 performed the best based on the cumulative ranking metric (Table (Table4).4). Conversely, the standard reference genes ACTB, GAPDH, and HPRT1 did not perform as well. These results confirm a high constancy of expression for a subset of the lung-NEPS genes using an alternative method in a secondary set of samples.
Using a publicly available microarray dataset (GSE10072, ), differential expression analysis was performed using SAM . Results from SAM were compared using the dataset normalized by MAS 5.0 alone, versus the same dataset normalized by MAS 5.0 with scale normalization using the lung-NEPS reference genes represented on the microarray. We observed differences in the total number of differentially regulated genes, depending on the normalization protocol used. When MAS 5.0 + NEPS normalization was used, 5502 genes were identified as up-regulated in cancer, whereas 4798 up-regulated genes were identified using MAS 5.0 alone. With respect to down-regulated genes, 2543 were identified using MAS 5.0 + NEPS, whereas 3325 were identified using MAS 5.0 alone (Figure (Figure3A).3A). According to the Canonical Pathway Analysis [Ingenuity Pathway Analysis (IPA)], we observe slight differences in both the number and the significance of identified pathways between the two sets of differentially normalized microarray data (Additional file 8). For example, while both datasets identify pathways such as mitochondrial dysfunction and protein ubiquitination, analysis of the dataset normalized by MAS 5.0 + NEPS identifies pathways known to be important in lung cancer, such as Neuregulin and JAK/Stat , at a higher significance relative to analysis of the same dataset normalized by MAS 5.0 alone (Figure (Figure3B).3B). Similarly, when evaluated using an additional publicly available lung cancer microarray dataset , we observe slight differences between the various pathways identified from analysis of differentially expressed genes derived from a NEPS-normalized dataset versus a dataset not normalized using the lung-NEPS genes (Additional file 9, Additional file 10). These results demonstrate that the choice of reference genes used for data normalization can influence the conclusions derived from gene expression studies.
In this study we present a methodology based upon permutation test analysis of SAGE data, to identify reference genes that more stringently satisfy the constancy requirements crucial for accurate normalization between samples utilized in gene expression experiments. Specifically, we have identified reference genes more effective for normalization than the traditional and previously reported housekeeping genes for lung, breast, and brain cancer gene expression profiling. Furthermore, we strongly emphasize that reference genes utilized for expression profiling should be selected in a tissue specific manner. Given that this methodology utilizes sequence-based data, its utility will increase as data generated from new next-generation sequencing platforms accumulate. The usage of more appropriate reference genes will have an impact on the interpretation of existing microarray data as well as expression data generated in future studies, and potentially will shed new insight into the molecular biology of cancer.
The authors declare that they have no competing interests.
RC analyzed the SAGE and quantitative PCR data to identify the key findings, and wrote the manuscript. KML led the construction of all lung SAGE libraries, contributed to interpretation of data and manuscript editing. LAP performed the quantitative PCR experiments and data interpretation. BPC and CEM provided insight into the statistical analysis as well as manuscript editing. CQZ contributed to the experimental work. THWC contributed to the data analysis. SL isolated the clinical samples, contributed to interpretation of results. CEM, MST, SL, RTN and WLL are the principal investigators on the project. They contributed to the overall design and data interpretation. All authors have read and approved the final manuscript.
The pre-publication history for this paper can be accessed here:
Summary of lung SAGE libraries. Raw tag count summary for lung SAGE libraries used in analysis.
Summary of breast SAGE libraries. Raw tag count summary for breast SAGE libraries used in analysis.
Summary of brain SAGE libraries. Raw tag count summary for brain SAGE libraries used in analysis.
Taqman probe IDs of genes assessed by quantitative RT-PCR. Applied Biosystems Taqman probe IDs for genes assessed by qRT-PCR.
Probes used for NEPS genes represented on the Affymetrix U133A microarray. Probes for NEPS genes represented on the Affymetrix U133A microarray.
Probes used for NEPS genes represented on the Agilent microarray. Probes used for NEPS genes represented on the Agilent microarray.
Permutation scores of previously identified and NEPS identified reference genes. Data used in scatter plot shown in Figure Figure11.
Ingenuity Pathway Analysis using genes from the analyses of a NEPS-normalized and unnormalized dataset by Landi et al. Ingenuity Pathway Analysis using genes from the analyses of a NEPS-normalized and unnormalized dataset by Landi et al.
SAM and pathway analysis of an Agilent lung cancer microarray dataset normalized with and without lung NEPS genes. SAM and pathway analysis of a dataset normalized with and without lung NEPS genes. (A) Number of probes identified as differentially over and underexpressed between cancer and normal using SAM on the dataset with and without NEPS normalization. Venn diagram illustrates the overlap in the genes identified as well as those which are different between the two analyses. (B) Canonical pathway analysis using Ingenuity Pathway Analysis. Dark blue bars represent the results from the dataset normalized with NEPS and median normalization and light blue bars represent the results from using median normalization alone. While similar pathways are statistically significant, each pathway is slightly different in the degree of statistical significance. Such differences illustrate the impact of reference gene selection and normalization on differential gene expression analysis.
Ingenuity Pathway Analysis using genes from the analyses of a NEPS-normalized and unnormalized dataset by Boelens et al. Ingenuity Pathway Analysis using genes from the analyses of a NEPS-normalized and unnormalized dataset by Boelens et al.
We thank Drs. William W. Lockwood and Ian M. Wilson for useful discussion and editing. This work was supported by funds from Canadian Institutes of Health Research. RC is supported by scholarships from the Michael Smith Foundation for Health Research and Canadian Institutes of Health Research.