Selection of endogenous control genes
We made a definition of established endogenous controls as genes available commercially, such as from Applied Biosystems. At the time of the analysis, they were: ACTB, B2 M, GAPDH, GUSB, HPRT1, PGK1, PPIA, RPLP0, TBP, and TFRC. From these, GAPDH, B2 M, PPIA, RPLP0, and TBP were selected as endogenous control candidates. They were selected, as described in methods.
The Ct value of each of these genes was submitted to the geNorm plugin for investigation of the stability index. The most stable gene pair was GAPDH and RPLP0. In order of decreasing stability, they were followed by TBP, PPIA, and B2 M. The exact definition of this method of classification is further described by Vandesompele et al. [2]. Briefly, the two most stable genes are identified by calculating expression ratios, over all samples, for all pairwise combinations of genes. For each pair of genes, the standard deviation over all samples is calculated, and for each gene, the mean of all standard deviations is calculated. Genes for which this value is highest are iteratively eliminated--in other words, the algorithm searches for genes that show the same expression profiles across all samples.
Comparison of real-time PCR expression data with microarray expression data
Real-time PCR data were obtained for 15 different genes of interest in addition to the 5 candidate endogenous controls. We compared the measurement of gene expression using these two methods of quantification. For each gene, this was done by creating a scatter plot, such as the example shown in Figure 1. Scatter plots for all combinations of genes, probe sets, and endogenous control combinations are found in Additional file 1. Pearson correlation coefficients were calculated for all of these combinations, and a summary can be seen in Figure 2. This figure prompts two discussions: a row-wise discussion of the differences in correlations between different genes and a column-wise discussion of the differences in correlations between different choices of endogenous controls.
Differences in correlations between genes
Some genes show good correlation between real-time PCR measurements and microarray measurements, and others do not. These differences can be explained biologically and technically.
It was investigated whether there was any systematic technical bias of the real-time PCR-to-microarray correlation. This was done by comparing the correlation metric to the mean absolute expression level and the standard deviation, both for microarray values and real-time PCR values. This is shown as scatter plots in Additional file 2, and no patterns could be identified.
Technical imprecision in individual measurements can also be a problem. All microarray scans were subjected to standard quality control measures, as detailed in the methods section, in order to exclude problematic samples. For real-time PCR, replicate measurements can clarify if deviance is a result of technical imprecision. A technical variation threshold is described in methods. For most genes, there were only a few measurements (< 10%) with technical variation above this threshold. The exception was ALOX12, which had 69 samples with coefficients of variance above this threshold. ALOX12 was present in very low quantities, with a mean Ct of 37.08 (Additional file 3). The discrepancy between measurements of this gene is therefore concluded to be due to technical imprecision in real-time PCR, as it is reflected in the lack of correlation between microarray and real-time PCR data. The data are included in the analysis to show the effects of technical imprecision, but its removal does not change the final conclusion, as specified in the sensitivity analysis in the methods section. In addition to the threshold values in the bottom of Additional file 4, we have included the raw data from the real-time PCR measurements as Additional file 3. This threshold could possibly be used as a measure of comparison for use in other experiments that do not have microarray data to compare with.
Because of alternative splicing and the difference in probe location for the two measurements, biological variability can also explain differences in measurement techniques for some genes. It is therefore not necessarily sufficient to talk about measuring a particular gene with two methods; the measurements should also be performed on the same part of the gene. The locations of real-time PCR probes and primers and of microarray probe sets are shown for the examples of IGF1 and EDNRA in Figure 3. For all other genes, the same plots are available in Additional file 4. ENDRA shows a high degree of correlation, even though different sections of the gene are measured with the two techniques - the opposite is true for IGF1. For IGF1, alternative splicing is unlikely to be the reason behind the lack of correlation, but unfortunately, for almost all other genes (11 out of 15), the real-time PCR primers and microarray probe locations are not in the same exon. Without specific real-time PCR measurements of relevant regions, it is therefore difficult to make a decisive conclusion on the gene-wise differences in correlation between microarray and real-time PCR. Focus shall therefore be on the column-wise discussion, as follows.
Differences in correlations between micorarray data and different real-time PCR normalizations
In Figure 2, the differences in correlations between different choices of endogenous controls can be read in each column. As demonstrated, the best correlation between microarray and real-time PCR is, on average, obtained when no endogenous control is used. The results obtained by using the geNorm method of Vandesompele et al. ranks 13th (for the GAPDH/RPLP0/TBP combination) and 17th (for GAPDH/RPLP0) among all 32 possible combinations of one or more of 5 endogenous control genes. This is an unexpected result, and it prompts a thorough discussion of the assumptions behind normalizations of PCR gene expression data.
A functional definition of gene expression should be expressed in units of mRNA per cell or nucleic mRNA concentration or a similar measure, because this is the level at which changes will affect the biology of the cell. In effect, the methods we compare here all rely on one of three assumptions to get to this functional definition: 1) that there exist one or more genes for which the expression is at a sufficiently constant level to be proportional to the amount of cells across all samples; this is the assumption behind the concept of endogenous control in ΔΔCt analysis; 2) that the overall distribution of the expression levels of all genes is the same in all cells, across all samples; this is the assumption behind the quantile normalization, which is the normalization part of the RMA pre-processing algorithm used for microarray analysis; and 3) that the total concentration of RNA is the same in all cells across all samples; this is the assumption behind the "no endogenous control" calculation of gene expression for real-time PCR data. This assumption stems from the fact that equal amounts of RNA were used in the sample preparation for the real-time PCR measurements, and it could as well be labeled "normalization to total RNA mass." It carries the further assumption that the reverse transcription of the RNA to cDNA does not introduce a bias.