Identification of a novel internal control through variation in gene expression levels

Differentially expressed genes, identified by microarray through global normalization, require validation of their expression patterns through Q-RT-PCR, which generally employees one internal control for normalization. This distinct normalization method calls for a correct internal control for the microarray and Q-RT-PCR data comparison. To prioritize the potential internal controls for Q-RT-PCR analyses and to study the possibility of general utilization of these potential internal controls, we first applied the block bootstrapping technique to rank genes with variance and IQR (Fig. 1). This method was then used to evaluate the gene expression patterns of the various internal controls in a lung adenocarcinoma microarray dataset (referred to as NHRI dataset) with the goal to provide insights into which internal controls might be a best choice in this study. For illustration, we first used the RMA for probe-cell level normalization and expression-value summary of the 66 microarrays in the NHRI dataset. Associated with each gene, the appropriate averages of 1,000 bootstrap replicates for the 5 computed statistics produced a box plot that disclosed variation in gene expression. As a result, many potential internal controls were revealed. To prioritize the targets, we removed those with lower intensity signals and the multiple probe sets on Affymetrix chip with potential cross hybridization contamination between genes as determined by NetAffx™ Analysis Center [20]. The remaining 14 candidate internal controls, which exhibited small variance expression, were identified using the box plots and these were DDX5, PKM2, BHLHB2, GLO1, LAPTM4A, SET, KIAA0152, CLTC, MSN, ABCF1, EPHB3, CCL5, PTPN21 and DDR1 (Fig. 2A). In addition, the box plots for 10 well-known internal controls [15, 21] (containing 23 probe sets) are also shown in Figure 2A. Surprisingly, 13 out of 23 probe sets have various levels of cross hybridization as identified by the NetAffx™ Analysis Center (Table 1). The unique characteristics of these well-known internal controls for Q-RT-PCR are their high expression levels and there are also huge variances.

Gene symbol	Probe set ID	Cross-hybridization	Column # in Fig. 2A
ACTB	AFFX-HSAC07/X00351_M_at	no	1
	AFFX-HSAC07/X00351_5_at	no	2
	AFFX-HSAC07/X00351_3_at	no	3
	200801 _x at	yes	4
GAPDH	AFFX-HUMCAPDH/M33197_M_at	no	5
	AFFX-HUMGAPDH/M33197_5_at	no	6
	AFFX-HUMGAPDH/M33197_3_at	no	7
	217398_x_at	yes	8
	213453_x at	yes	9
	212581_x_at	yes	10
B2M	216231_s_at	yes	11
	201891_s_at	yes	12
HMBS	203040_s_at	no	13
HPRT1	202854_at	no	14
RPL13A	212790_x_at	yes	15
	211942_x_at	yes	16
	210646_x_at	yes	17
	200716_x_at	yes	18
	200715_x at	yes	19
RPL32	200674_s_at	no	20
SDHA	201093_x_at	yes	21
UBC	211296_x_at	yes	22
YWHAZ	214848_at	no	23

The copy number of the individual housekeeping gene chosen for relative quantification should be in a similar range to that of the target gene to make comparative quantification possible [22]. Further analysis indicates that we can identify a series of internal control candidates, which have characteristics of small variance in different microarray intensity intervals (Additional file 1). These potential internal controls are presented in different intensity ranges in order to appropriately normalize different target genes. Despite the fact that these potential internal controls may exhibit a greater variance than DDX5 or other internal controls listed in Figure 2A, these potential internal controls have much smaller variance than ACTB and GAPDH (Additional file 1, right portion). This finding further supports a view that the intensities of normalized microarray data and the copy numbers of Q-RT-PCR detections in gene expression patterns could be examined in a similar range. Our approach may provide a method to identify potential internal controls to be in a similar range of expression as the selected target genes.

To prioritize the potential internal control for the lung cancer microarray data, two major public accessible lung cancer microarray datasets, which also used Affymetrix chips, namely the Boston and the Ann Arbor datasets [23, 24], were included for data comparison. Eleven candidates, after the exclusion of ABCF1, BHLHB2 and LAPTM4A, were available in both Boston and Ann Arbor datasets and were included in the analysis. Figure 2B and 2C show the results of basic bootstrapping using these two datasets of unpaired design. All 11 candidates exhibited less variation than most well known internal controls, suggesting that all 11 candidates have potential to serve as an internal control, at least for lung cancer. To finalize the target for further empirical validation, these potential controls were sorted in the order of increasing gene expression intensities and decreasing IQR, respectively. As a result, CLTC, DDX5 and MSN were found to exhibit sufficient intensity. However, DDX5 gave the smallest variation among the three. Therefore, in this study, we chose DDX5 for further characterization.

Considering this study was a paired design with a moderate sample size, we applied the block bootstrapping technique and evaluated the gene expression patterns of various internal controls using the NHRI dataset (as described in Materials and Methods). A total of 39 blocks, resulting from 27 pairs of adjacent normal and tumor samples and 12 un-paired samples of microarrays, were used in the block bootstrap. The re-sampling process was repeated 1,000 times to obtained 1,000 bootstrap replicates of the minimum, first quartile, median, third quartile and maximum for the expression level of each gene. The bootstrap replicates for all five statistics for DDX5 expression are roughly constants, but those for GAPDH expression vary greatly (Fig. 3).

DDX5 exhibited relatively similar expression patterns between adjacent normal and tumor samples as well as across various lung cancer cell lines when compared to GAPDH (Fig. 4). In addition to RMA, we also employed the other two different methods, MAS5 and GC-RMA, for the probe normalization to analyze 66 Affymetrix microarray chips using lung adenocarcinoma RNA samples. The three different normalization methods displayed consistent expression patterns for DDX5, suggesting that DDX5 may serve as a reliable internal control [25].

Q-RT-PCR validation and comparison with microarray

The expression pattern comparisons of DDX5 in other microarray datasets

To explore the possibility that DDX5 can be generally utilized as an internal control, we further analyzed DDX5 expression patterns using a new Affymetrix based NCI60 dataset [27] from Genomics Institute of the Novartis Research Foundation (GNF). The expression patterns of DDX5 were investigated in 8 different types of cancer cell lines (NCI60) and several additional cell lines with a total of 84 cell lines and DDX5 exhibited much smaller variation than GAPDH for all cell types (Fig. 5A). Moreover, the expression patterns of DDX5, CLTC and MSN (Fig. 5B) were further compared with ACTB and GAPDH in the lung cancer cell lines found in NCI60. All three candidates showed less variance than both well-known internal controls. Similar patterns were also observed in our NHRI dataset (Fig. 2A). Additional cDNA microarray datasets were also downloaded from Stanford Microarray Database [28], including lung cancers [29], hepatocellular carcinomas (HCC) [30] and the HeLa cell cycle dataset [31], which reported cell cycle related genes by synchronizing HeLa cells at different phases of cell cycle [32] (Fig. 5C). Based on these comparisons, we conclude that DDX5 is evidently a novel internal control with a relatively constant expression pattern in lung adenocarcinoma. In addition, using limited microarray datasets analyzed in this paper, the variation in DDX5 levels was also relatively small in other cancer cell lines, raising the possibility that DDX5 could serve as a novel internal control.

Sample preparation

A total of 66 samples were used for microarray analysis, including paired adjacent normal-tumor samples from 27 patients underwent surgery for lung cancer at the Taipei Veterans General Hospital, two tissue mixtures from the Taichung Veterans General Hospital (one was adjacent normal lung mixtures and the other was lung adenocarcinoma mixtures), two commercial human normal lung tissues (Clontech (Catalog No. 636524) and Stratagene (Catalog No. 735020)), one immortalized, nontumorigenic human bronchial epithelial cell line (NL-20 (ATCC^® No. CRL-2503™)) and 7 lung cancer cell lines (A-549 (ATCC^® No. CCL-185™), NCI-H1299 (ATCC^® No. CRL-5803™), NCI-H661 (ATCC^® No. HTB-183™), CL_1-0 [33], CL_1-1 [33], CL_1–5 [33], and CL_1–5-F4 [34]).

RNA extraction and reverse transcription

We used the total RNA samples for Q-RT-PCR analyses. RNA preparation and analysis were performed according to the previous study [11]. Briefly, the quality of the total RNA for microarray analysis was determined using Spectra Max Plus (Molecular Devices) and had an OD260/OD280 ratio ranging from 1.9 to 2.1. RNA was subjected to reverse transcription with random hexamer primers. To hydrolyze contaminating DNA in the RNA preparations, RNA was incubated with amplification-grade DNase I (Life Technologies, Gaithersburg, MD). After incubating the reaction mixture, the reaction was stopped by heating at 65°C. After DNase treatment, the RNA was subjected to reverse transcription reaction by the ThermoScript™ RT-PCR system (Life Technologies) and cDNAs were then used in the Q-RT-PCR.

Q-RT-PCR

Microarray experiments, normalization and Wilcoxon signed-rank test

Protocols, reagents for hybridization, washing and staining followed previous methods [35] and the Affymetrix's instructions [36]. Labeled cDNA was hybridized to the Affymetrix GeneChip Test 3 Array to verify quality prior to hybridize to the Affymetrix Human Genome U133A Array. The data discussed in this publication have been deposited in NCBIs Gene Expression Omnibus (GEO) [37] and are accessible through GEO Series accession number GSE7670. The images were transformed into text files as intensity information using the MAS5.0 (Microarray Suite software 5.0) developed by Affymetrix [38, 39]. We used three array processing methods to produce and normalize the Affymetrix expression signals for the transcripts based on corresponding probe pairs of oligonucleotides. MAS5.0 as provided by Affymetrix was used to carry out the probe-pairs adjustment. Both MAS5.0 and RMA [40], Robust Multichip Average via quantile adjustment, are commonly used. Another algorithm GC-RMA [41] is similar to RMA, but pools the probes with comparable numbers of G-C bonds to achieve a stable mismatch adjustment. Based on paired adjacent normal and tumor samples of 27 patients, we used the Wilcoxon signed-rank test [42] to identify differentially expressed genes. In contrast to Student's t-test inappropriately used in some studies [43], the signed-rank test is distribution-free and adjusts for the paired design.

Block bootstrapping

To assess the variation in the microarray expression level for a specified gene in an experiment of moderate sample size and possibly including paired samples, we designed a block bootstrapping procedure to analyze the microarray data from 66 lung samples, containing 27 pairs of patient adjacent normal-tumor samples and 12 un-paired samples. Bootstrapping is a breakthrough statistical approach using a computationally intensive re-sampling technique and it allows complex problems to be solved in which the accuracy of a devised statistical procedure can not be analytically evaluated [44, 45]. Block bootstrapping was originally named for re-sampling methods in dependent cases, especially time series data [46]. The basic bootstrap generates artificial samples that allow the making of an inference of interest through re-sampling the original data with replacement in which all observations are assumed to be mutually independent and from the same distribution. To guarantee the structure of independence in bootstrap re-sampling, we employed the concept of blocking for the paired data by treating each individual patient, mixture tissue or cell line as a block. By selecting an observation within the block with equal probability when combined with all the other un-paired samples, we obtained an independently re-sampled dataset. We then created a bootstrap sample by randomly sampling the blocks in the dataset, and computed the bootstrap replicates of the relevant summary statistics of the expression levels. Repeating this bootstrap re-sampling scheme sufficient times, such as 1,000, we then used the averages of these bootstrap replicates to reveal the variation in expression summaries corresponding to a specific gene across the microarrays. Appropriate internal controls can be selected by ranking the variations in gene expression.

Correlations of microarray and Q-RT-PCR data

To explicate the correlation between microarray and Q-RT-PCR in this study of paired design, we calculated the differences between the log-scaled measurements of the Q-RT-PCR and microarraydata from the tumor and adjacent normal tissues of 24 patients. The other 3 paired samples did not have sufficient materials for further Q-RT-PCR analysis. Pearson's and Kendall's τ correlation coefficients were then tested at a significant level of 0.05. Summarization of the expression levels and normalization for microarray data were conducted using GeneSpring^® 7.3 (Silicon Genetics, Redwood City, CA). The computer programs for the block bootstrapping method and correlation using R 2.1.1 [47] are presented in the Additional File 2.