Selection and validation of endogenous reference genes using a high throughput approach

Data processing

The GenePix Pro 4.0 software was used for array image analysis and calibration. All statistical analyses were performed with the SPLUS package http://www.insightful.com. Thirteen arrays with missing data in more than 30% of the spots and 8 arrays with irregular distribution according to M-A plots (M and A representing respectively log-ratios and average intensities) [32] were excluded. Spots in which > 50% of the pixels reached saturation in either channel, flagged spots and spots with intensity ≤ 200 in one channel and ≤ 500 in the other channel were filtered out. If the intensity in one channel was < 200 but that of the other channel was > 500, we arbitrarily assigned an intensity of 200 to the channel with the lowest intensity. The _log2Cy5/Cy3 ratios were normalized by approximating median values to zero. Spots with ratios > 100 or < 1/100 were truncated at 100 and 1/100 respectively. Data normalization was done by median centering _Log2ratio of all genes for each array. For data normalization we used certain arbitrary cutoff criteria to remove spots with weak signal. Weak signal approximating background fluorescence is not reliable and the corresponding log ratios are often distorted resulting in disproportionately high or low values that would bias the statistical results. Because the interpretation of spots with low signal is difficult to make we have adopted the policy of excluding them to focus the analysis on genes whose expression pattern is more reliably tracked by the array tool. In regard to spots with a ratio of >100 or <-1/100, most of them are due to extremely weak signal in one channel which generates a disproportionately extreme value. Although there is no published "gold standard" for the selection of such cutoffs, in practice the parameters that we used for this study are most commonly accepted by investigators and reasonable because allow retention of most of the data in the array excluding only the most extreme and least reliable information.

Validation of the data set through analysis of "reference concordance"

To select from the data base the most reliable genes, we first assessed reproducibility by testing the level of concordance in 14 reference experiments. Reference concordance relies on the expectation that results obtained by repeatedly hybridizing the same test and reference material should perfectly collimate and the degree of deviation from such prediction estimates experimental variance. Concordance can be easily measured by periodically re-hybridizing the same test sample with the constant reference sample. We analyzed a matrix of 7 forward and 7 reciprocally-labeled replicate array experiments hybridized periodically every other 25 cDNA array slides. Reciprocal labeling was applied to measure labeling bias [30]. Data generated by reciprocal experiments were mathematically reversed into the same labeling direction for data analysis. Genes that were discordant due to labeling bias were identified by student's t test as those with ratios significantly different between the 7 forward and 7 reciprocal experiments after reversal of the reciprocal values (P < 0.05) and the median ratio difference between the 7 forward and 7 reciprocal experiments was larger than 1.5 fold. The genes whose variances of the log ratio across all experiments were among the top 1 percentile of variance of all genes were identified as discordant due to hybridization bias. In addition, genes with more than 50% missing values were excluded. Overall 1,343 out of 16,738 were judged potentially discordant and were excluded from further analysis. The remaining 15,395 concordant genes were used for subsequent analyses.

Endogenous reference gene identification

We studied 384 of 419 consecutive array experiments remaining after the exclusion of the 14 reference concordance experiments and 21 experiments judged of poor quality due to missing data or irregular distribution by M-A plot analysis. Selection of candidate endogenous reference genes followed these steps.

Among these genes, we further selected candidate endogenous reference genes according to the following criteria: in at least 95% of the 384 arrays the _Log2 ratio had to be within 2 SD from zero and the average intensities of both channels across all samples need to be higher than 2,000. The intensity parameter was added to ensure that the selected endogenous reference genes were expressed at a relatively high level in most tissues to mitigate excessive fluctuations in _Log2 ratios occurring when a low value is applied as denominator. This is important when a reference gene is applied as a denominator in the equation used to normalize the ratio of other tests genes; in such cases robust gene expression decreases the range of ratios resulting from the analysis decreasing, therefore, the experimental variance. Sixty-nine genes fit these criteria. The range of mean _Log2 ratio was from 0 to 0.15 and the SD of mean _Log2 ratio was from 0.27 to 0.42. Analysis of 7 replicate array based on the same group of candidate endogenous reference genes demonstrated a similar distribution of mean _Log2 ratio and SD of mean _Log2 ratios (with exception of TFAP2B) suggesting that this variance could be attributed to predominantly experimental noise. Among these genes we further selected for validation 11 that had high mean fluorescence intensity (underlined in Table 1). However, probe preparation and other technical considerations limited to analysis to only 8 of these genes which included: NEED8, HIRIP5, GPLD1, RAB31, SNRPD3, FZD6, COL8A1 and AFAP.

"Leave-one-out cross-validations" were used to estimate the error [33, 34]. In each validation experiment one array was left out. The remaining 383 arrays were used to identify endogenous reference genes using the same parameters used for the original analysis (test/reference ratios < 2 SD in at least 95% of the experiments and average intensity > 2,000 in both channels). The median _Log2 ratio of the endogenous reference gene in the left out array was then compared with the "real" normalization factor consisting of the median _Log2 ratio for all genes in that array supposed to approximate a balanced expression of test and reference genes for that array. This procedure was repeated for each array 384 times. Only in 3 arrays out of 384 (2.3%) were found errors in endogenous reference gene selection (errors were considered median _Log2 ratio of endogenous reference genes in the left out array > 2 standard deviations apart from the "real" normalization factor in the same array (which is zero for normalized arrays).

Test gene expression estimates according to house-keeping gene selection: the transcriptome array

To validate the usefulness of candidate endogenous reference genes in large sample populations we developed a new tool that we are planning to use in the future for validation of gene expression across extensive data sets. This tool displays cDNA libraries originated from different tissues or cell lines individually spotted on a solid surface. The principle of this technology is similar to RNA dot blot which uses RNA isolated from samples and transferred to membrane making, therefore, the transcriptome array a high-density dot blot. The labeled DNA probe of interest is hybridized to the immobilized complementary strain of mRNA. A reference gene hybridization will carried out simultaneously to estimate the relative expression of the gene of interest compared with the reference gene. The new tool we describe here, termed transcriptome array, utilizes cDNA generated from source mRNA for target immobilization to improve the stability of the immobilized targets and differentially fluorescence-labeled test and reference probes (RNA or double strained DNA) then can co-hybridize on to the same spots. Using a validated RNA amplification technology [30], large quantity of pure amplified RNA with proportional representation of source mRNA species could be generated from which cDNA could be obtained through a reverse transcription reaction. Because of the minimum amount of cDNA used for fabricating each transcriptome array (<5 nano gram cDNA/spot) and the size of spots (100 um), the expression of a large number of genes can be theoretically analyzed on thousands of different samples simultaneously. Since the amount of cDNA spotted may vary according to the quality of the starting material and the efficiency of RNA preparation for each sample, absolute estimation of fluorescence from the hybridized probe is not informative of the expression of the given gene in each sample. Therefore, a reference system is necessary so that comparative expression of the test gene can be presented proportional to that of a consistently expressed gene across all samples. Therefore, interpretation of data derived from the transcriptome array relies on endogenous reference gene normalization.

To test the usefulness of various endogenous reference gene we resorted to the well characterized expression of the melanocytic lineage-specific gene gp100/Pmel17 [31] that is expressed exclusively though heterogeneously in cells of melanocytic lineage [35, 36]. The assumption of this experiment was that in spite of its heterogeneity of expression in samples of melanocytic origin, overall the expression of gp100 should be higher in meloanocytic compared with non melanocytic samples with a high degree of significance. Rho-C has been associated with the metastatization process in melanoma but its specificity for melanocytic lineage is unknown [37]. The differentiation control element DICE is found in the 3'-UTR of numerous eukaryotic mRNAs and there is no solid association between its pattern of expression and specific physio-pathological or developmental conditions [38]. Therefore, we compared the expression profile of these three genes in 829 cDNA libraries that included 106 melanoma cell lines, 127 melanoma metastases, 2 benign nevi (total of 235 melanocytic samples) and 593 miscellaneous samples containing a predominance of primary esophageal, renal and colon cancers paired with normal tissues from the same organ and a large number of circulating mononuclear cells. In addition, samples from most other normal and cancerous tissues were included although in smaller number (complete list of samples available upon request).

As endogenous reference genes we chose β-Actin and new candidate genes identified by this study (see previous section). The test genes (gp100, Rho-C and DICE) were separately hybridized to the transcriptome array. Each gene labeled with Cy3 was co-hybridized with individual endogenous reference genes labeled with Cy5. We then compared _Log2 Cy5 / Cy3 ratios between melanocytic and non-melanocytic tissues. gp100 was, as expected, expressed more in melanocytic lesions with high degree of significance no matter what gp100 / endogenous reference gene combination was used. The range of significance, however, varied extensively depending upon the gp100 / endogenous reference gene combination. This difference was considered representative of the ability of different endogenous reference genes to normalize for tissue-specific gene expression patterns (Figure 1). In details, all endogenous reference genes could segregate melanocytic from non melanocytic lesions with a high degree of significance (Figure 2; unpaired two sample t test p₂-values ranged from 6 × 10^-8 to 3 × 10^-36). There was, however, a big range in the discriminatory capacity among endogenous reference genes with NEDD8, RAB3 and FZD6 providing the highest predictive value (t test p₂-value = 3 × 10^-36, 1 × 10^-32 and 1 × 10^-21 respectively) and β-actin being one of the least reliable (t test p₂-value = 1 × 10^-8).

The large variation in the results obtained using different endogenous reference genes to normalize gp100 expression could have been due to a higher stability of the expression of some genes across all samples or to a differential expression of the endogenous reference genes themselves in melanocytic lesion. In the latter case, the better results obtained could be coincidental and not useful in other experimental situations. We, therefore, tested whether pooling results obtained with all endogenous reference genes independently of each predictive value in this controlled experimental situation could yield results as informative about gp100 pattern of expression as those obtainable with individual endogenous reference genes, particularly those that provided the best prediction. The advantage of this strategy is that it does not depend on previous knowledge of gene expression patterns for the selection of individual endogenous reference genes not applicable in conditions, in which the suitability of a gene for a given experimental situation is, contrary to gp100, is unknown. Thus, we compared the predictive value of the average of gp100 / endogenous reference gene _Log2 ratios obtained with all endogenous reference genes and those obtained using the genes that provided the best five or the best three results (Figure 2). The use of pooled information from various endogenous reference genes appeared to stabilize the non-melanocytic _Log2 ratios gp100 / endogenous reference. In fact, the SD of the_Log2 ratios gp100 / endogenous reference among non-melanocytic samples derived by pooling endogenous reference gene results were significantly lower than the SD obtainable with any individual endogenous reference gene (F-test). The basis for this test was that non-melanocytic lesions uniformly should not express gp100 and, therefore, the _Log2 ratios for any gp100 / endogenous reference gene combination should be constant resulting in very low SD. SD for non-melanocytic lesions were much lower using pooled endogenous reference genes for the normalization (Figure 2) and this was true at extreme levels of significance (F-test value for SD when all endogenous reference genes were used for normalization compared with individual endogenous reference were 0, 1.5 × 10^-9, 0, 5.6 × 10^-7, 9.0 × 10^-5, 0, 0, 0 and 0 for NEDD8, RAB3, FZD6, AFAP, COL8A, SNRPD, HIRIP5, β-actin and GPLD1 respectively). Similar results were obtained when results from the best 5 and best 3 endogenous reference genes were pooled together. In addition, even if pooling data together did not provide the highest level of discrimination compared to the best results obtained with selected individual endogenous reference genes, the capacity of pooled data sets to discriminate gp100 expression between melanocytic and non melanocytic samples was still very high (2 × 10^-22; 4 × 10^-29 and 1 × 10^-34 respectively when all, the best 5 or the best 3 endogenous reference gene results were pooled together). Interestingly, different endogenous reference genes provided not only different levels of discrimination but also provided different estimates of gp100 expression in melanocytic lesions with gp100 / endogenous reference gene _Log2 ratio above (NEED8, RAB3, FZD6, AFAP, SNRPD3, HIRIP5, GPLD1) or below (COL8A1, β-actin) 0 (Bottom graph, Figure 2). This, of course, has no impact on the ability to characterize patterns of expression of different genes in different tissues but rather suggests that the latter group of endogenous reference genes is expressed at higher copy number, therefore, biasing the data toward its own fluorescence channel. Simple adjustment in probe concentration could easily solve this problem.

We then analyzed the expression of Rho-C and DICE whose pattern of expression in the transcriptome array cannot be predicted based on available information save for the notion that of the expression of Rho-C is related to the metastatic process of melanoma [37]. The pattern of Rho-C was characterized with all the endogenous reference genes and representative information is presented in Figure 3. Most endogenous reference genes yielded a pattern suggestive of a preferential expression of Rho-C in melanoma metastases. In addition, averages of Rho-C / endogenous reference gene _Log2ratios demonstrated clearly a specific expression of this gene in melanoma metastases. This pattern was not necessarily specific for melanoma as the other cancerous tissues spotted in the transcriptome array were obtained from primary lesions. Thus, it is possible that the preferential expression in melanoma metastases was due to the metastatic process rather than the melanocytic lineage. Interestingly, the melanoma cell line Mel-A375 did not constitutively express Rho-C as previously observed by Clark et al. [37]. In addition, none of the melanoma cell lines mostly derived from melanoma metastases demonstrated expression of Rho-C. This information suggests that Rho-C may be involved in the natural metastatic process in vivo but is not constitutively expressed in vitro. Finally, DICE did not demonstrate any specificity of expression and appeared variably expressed in all specimens independently of the endogenous reference gene used for normalization.

Validation of data by protein analysis and quantitative real-time PCR (qRT-PCR)

We tested whether gp100 / endogenous reference gene _Log2 ratios correlated with gp100 protein expression level as measured by intra-cellular FACS analysis (Figure 4). The gp100 / endogenous reference gene _Log2 ratios were derived from the transcriptome arrays while protein expression was based on previous characterizations of the same cell lines [35, 39]. Overall, a good correlation (Pearson's correlation) was noted for gp100 values normalized with most endogenous reference genes. However, with the exception of NEDD8, the best correlation was obtained when pooled information was used from the best 5 and 3 candidate endogenous reference genes as defined before. Interestingly, RAB3 that scored very high as a predictor of gp100 transcriptional expression in melanocytic lesions and β-actin provided the worst correlation with protein expression. This data suggest that, although different endogenous reference genes may yield better predictive value than others in large data sets, it is likely that their reliability varies in different conditions and possibly the best results can be obtained by pooling information from several of them.

We also tested by quantitative real-time PCR (qRT-PCR) the discriminatory capacity of endogenous reference genes in samples with known patterns of expression of gp100. We selected 6 previously characterized [35, 36] melanoma cell lines four of which were known to express (MEL-553B, MEL-1317, MEL-526, MEL-1102) and two known not to express (MEL-836, MEL-1195) gp100 at the transcriptional and protein level (Figure 5) [36, 39]. In addition, 4 fine needle aspirate (FNA) samples from melanoma metastases expressing gp100 [40] and a series of tumors (2 renal cell, 2 esophageal one gastric carcinomas and HELA cells line) known not to express gp100 were tested. Absolute expression of gp100 was measured and although the lesions supposed to be gp100 expressing appeared more positive than those supposed to be gp100 negative there was no absolute demarcation among the two groups (Figure 5a). A better differentiation between positive and negative lesions could be obtained when the various endogenous reference genes were used for the analysis or the combined values from all of them (Figure 5b).

Target preparation for cDNA microarray hybridization

Total RNA extracted from test and reference samples was transcribed in vitro into anti-sense RNA (aRNA) and reverse-transcribed into fluorescence labeled cDNA for hybridization to 17,000 gene cDNA-based microarrays [30, 44]. Pooled PBMC were used to prepare reference aRNA to be co-hybridized in all experiments with test aRNA. cDNA targets were labeled with Cy3 (green) for reference material and Cy5 (red) for test material with the exception of reciprocal experiments.

cDNA Microarrays

cDNA (UniGene cluster) microarrays were printed at the Immunogenetics Section, DTM, CC, NIH with a configuration of 32 × 24 × 23 containing 17,500 elements. Clones used for printing were selected from the Research Genetics RG_HsKG_031901 8 k clone set and the 9k RG_Hs_seq_ver_070700 40 k clone set. The 17,500 spots included 12,072 uniquely named genes, 875 duplicate genes and about 4,000 expression sequence tags. For a complete list of genes included in the Hs-CCDTM-17.5k-1px printing please visit our web site at http://nciarray.nci.nih.gov/gal_files/index.shtml.

Transcriptome array

A collection of aRNA-based libraries was prepared from 960 frozen tissue samples or cell lines and individual aRNAs were reverse transcribed into cDNA in the presence of 1 μl of dN6 primer (8 μg/μl), 6 μl of first strand buffer, 3 μl of 10 mM dNTP, 3 ul of 50 mM dTT, 1.5 μl of Rnasin and 10 μl of aRNA (6–12 ug) with addition of 3 μl of Superscript II (BRL) in 5.5 μl volume after heating to 65°C for 5 minutes. Reactions were carried out in 96-well plates at 42°C for 90 minutes followed with addition of 1 μl Rnase H and heating for 20 min at 37°C. cDNA were further purified by utilizing CentriCept gelfitration plates to remove non incorporated primer and dNTP. Purified cDNA were transferred to 384-well plates and dried by speedvac. Samples were re-suspended in 13 μl of 3 × SSC followed by shaking at 3,000 rpm for 15 min. Reconstituted cDNA libraries from individual samples and spiked reference gene at different concentration were printed on to poly-L-lysine coated glass slides at the concentration of 0.5–1 μg/μl using the OmniGrid (GeneMachine. San Carlos, CA) printer and Telechem printing pins (TeleChem International, Inc. Sunnyvale, CA). Each 100 μm diameter spot was duplicated at a 250 μm distance. After complete exsiccation, slides were post-processed as described at http://nciarray.nci.nih.gov/reference/Post_Processing.html.

Probe design Transcriptome array hybridization experiments

Specific PCR primers for the amplification of each gene probe were designed using Primer 3 program http://molbio.info.nih.gov/molbio/analysis.htm. In order to generate single strand fluorescence labeled cDNA for hybridization, modified specific 5' primers with an extension of the T7 promoter region were used for PCR amplification. Double strand PCR fragments (300–800 bp) were used as template for in vitro transcription to generate single stranded mRNA. 3 μg of amplified RNA were then reverse transcribed into cDNA in present of 4 μl of first strand buffer, 1 μl dN6 primer (8μg/μl; Boehringer Mannheim), 2 μl 10 X low T-dNTP (5 mM A, C and GTP, 2 mM dTTP), 2 μl Cy-dUTP (1 mM, Cy3 or Cy5), 2 μl 0.1 M DTT, 1 μl RNasin, 3 μg amplified mRNA in 8 μl DEPC H₂O. After heating to 65°C for 5 minutes and cooling to 42°C, 2 μl of SSII was added and the labeling reaction was carried out at 42°C for 90 min. Probe purification and hybridization were performed as previously described [30, 44].

Quantitative real-time polymerase chain reaction (RT-PCR)