Background

The generation of very large amounts of gene microarray data poses a challenge, not only for its processing but also for its interpretation, due to intrinsic false discovery rates. In addition, the problem of background noise along with differences in hybridization efficiencies is also an important factor generating variability within and among microarray chips, constituting major confounding elements in gene expression analysis. Some of the most advanced commercially available software can automatically account for most, but not all, of these challenges.

In general, the persistence of confounding elements generates the need for appropriate data normalization methods, such as using the specific nth percentile signal intensity value of a particular array. Often, software-generated normalization methods have already incorporated the nth-percentile approach (e.g., 75th percentile) along with some background-subtraction mechanism [1]. In this regard, dealing with the assay’s inherent background noise becomes critical to account for signal stringency. Hence, the importance of using filtering parameters to accomplish the desired signal-to-noise ratios becomes obvious.

Moreover, it becomes necessary to use housekeeping or reference genes as endogenous controls for further gene signal normalization. This is not a common use in gene microarray analysis, where log transformation, background subtraction, and nth percentile normalizations have been the norm [1]. In this regard, the use of endogenous control genes is prerequisite in qRT-PCR and relative RT-PCR [2–4]. The principle behind this methodology consists of simply using widely expressed genes that do not respond to most treatments as references to compare to genes of interest (target genes) that do change. This helps in the proper interpretation of gene expression patterns and in calculating relative gene expression fold changes between treatment groups, minimizing technique-derived experimental errors. Thus, the same rationale should apply in gene microarray analysis. However, selecting the right endogenous control genes for normalizing data can be difficult since these widely expressed reference genes are not truly universal. In this regard, there have been observed changes in reference gene expression with different treatments as well as tissue-specific differential reference gene expression patterns [2–10]. For this reason, a systematic method must be determined for its use in the selection of array-specific (i.e., tissue- and taxa-specific) endogenous control genes based on a pool or pools of pre-established and widely used housekeeping or reference genes.

In the present examination, dependent measure data sets were derived from paired DNA microarray gene expression analyses performed on whole-blood samples from homozygous ZDF rats exhibiting clinically-relevant type 2 symptomatology in comparison to heterozygous healthy lean controls [11]. In this regard, the ZDF rat has been well established in the biomedical literature as a high-resolution translational model for elucidation of underlying pathophysiological mechanisms critically linked to advanced therapeutic development for major human disorders, including type 2 diabetes [12,13], cardiovascular disease [14], renal disease [15], atherosclerosis [16–18], and rheumatoid arthritis [11]. In addition, a list of potentially suitable endogenous control genes for the study of whole-blood ZDF rat samples is provided.

Material and Methods

STATISTICAL ANALYSES:

Software-generated data was compared using moderated t test method with Benjamini-Hochberg multiple testing correction. The t test was used for evaluating the manually normalized data.

Results

APPLICATION OF FILTERING PARAMETERS:

The first step before data analysis deals with signal quality control and the setting of filtering parameters in a particular data set. In this evaluation, the filters previously described were used to identify potential endogenous control gene candidates from the list of 34 widely used reference genes [2–8,19] (Table 1). After this initial filtering process, in which genes not meeting the signal quality criteria were filtered out (i.e., flagged as “compromised”; S/N <2), a working list of 18 gene probes corresponding to 16 endogenous gene candidates was made (Table 2). It should be noted that there can be more than 1 probe per gene, each having a different sequence, thus hybridizing to a different region of the gene transcript.

SOFTWARE-GENERATED GENE SELECTION: Suitable endogenous control genes should exhibit minimal-to-no expression variation between groups (e.g., control vs. treatment), in this particular study, between diabetic obese and healthy lean groups. In this evaluation, an absolute fold change value of 1.2 was set as the limit for the gene selection criterion. In this way, further filtering by fold change yielded 3 suitable endogenous control gene candidates, which can be seen in Figure 1. In addition, Table 3 shows these software-selected genes (Hsp90ab1, Pum1, and Srsf4) along with their fold change and p-values.

MANUAL NORMALIZATION: The manual method involving 75th percentile normalization of background-subtracted signals along with fold change calculations yielded 5 potentially suitable endogenous control candidates. Three were the same as the software-generated genes, plus 2 additional genes – Dimt1 and Gusb (Table 4). These genes exhibited fold change values <1.2 and >−1.2, with p-values considered not significant (p>0.05). In this regard, after manual normalization, there was an additional gene, Decr1, that exhibited an acceptable fold change value of 1.19 but had a p-value <0.05 and hence was not selected (Table 2). Table 5 shows an arbitrary gene grouping based on signal intensity values (low, medium, high). This helps in the selection of suitable endogenous control genes because, as mentioned earlier, these genes should ideally be chosen so that they encompass a large signal intensity spectrum in such a way that it compensates for the potentially diverse copy numbers of target genes (translated as signal intensities) [3]. Finally, Table 6 shows the Agilent probe sequences of each of the endogenous control genes selected by this method.

Discussion

FILTERING PARAMETERS TO ACCOMPLISH THE DESIRED SIGNAL-TO-NOISE RATIOS:

Gene microarray data need to be adequately filtered. The first step before data analysis involved a quality control step in which irregular signals or “compromised” features are removed (i.e., filter by flags: detected; not detected; compromised). Often, in order to accomplish an acceptable microarray signal intensity level, a signal should be at least twice as strong as that of the background (i.e., signal-to-noise ratio ≥2) and, depending on the desired stringency level, this filter cut-off can be set to a signal-to-noise ratio of 3 or higher. Normally, gene microarray technologies produce a consistent background signal whose mean level information can be easily obtained from the raw image data (e.g., using feature extraction software). A microarray platform’s software automatically subtracts calculated background signal from raw signal values, effectively yielding processed raw signal values. To filter out genes whose processed raw signal values are less than twice the background (S/N <2), a filter (i.e., processed raw signal cut-off) should be set to a lower limit equivalent to the array’s mean background signal value. In this way, if the processed (i.e., background-subtracted) raw signal is added to the technology’s mean background signal value, it will be equivalent to S/N=2. Processed data falling below the S/N=2 level were eliminated from the final data set. Moreover, this filtering by expression level was applied so it would accept a gene when in at least 1 of the subject groups studied (e.g., control; treatment “A”; and treatment “B”) is detectable since, for example, a given treatment/s could cause downregulation of a gene below a level corresponding to S/N=2. The same applies when a gene is only detectable after a treatment. In this way, when a particular gene or group of genes was present at S/N ≥2 in at least 1 of the subject groups, then those genes passed the filter.

ENDOGENOUS CONTROL GENE SELECTION AND NORMALIZATION:

As mentioned earlier, selecting the right endogenous control genes for data normalization can be difficult due to gene expression changing with treatments or to tissue-specific differential gene expression patterns [4,9,10]. For this purpose, based on several important publicly available studies [2–8,19], a list of 34 widely used reference genes was built to be evaluated with the experimental data.

Ideally, it is preferable to select more than 1 endogenous control gene and to average their values. In this regard, it is recommended that, when possible, the endogenous control genes be chosen so that they will exhibit different signal intensity levels (e.g., low, medium, and high) [3]. This would account for the differences in copy number (i.e., signal intensities) among the target genes. Hence, the signal intensity levels of the potential endogenous control gene candidates were compared to be sure they spanned a relatively wide range. Moreover, suitable endogenous control genes should be selected such that each is involved in a different cellular function and/or is found in different chromosomes [5,6]. Although this may not always be possible, it is recommended that at least 2 of these criteria be satisfied (Table 5).

Finally, in order to overcome or minimize inter-array differences, scaling to the nth percentile is recommended (in this case, to the 75th percentile) [1]. If using a linear scale (i.e., not log-normalized), as in this case, the processed (background-subtracted) signal intensity values are divided by the 75th percentile value corresponding to the particular array. This scaling is applied to both potentially suitable endogenous control genes and target genes.

NORMALIZATION OF TARGET GENES:

After finding and normalizing suitable endogenous control genes, the next step is to use them to normalize genes of interest to calculate their expression pattern though fold change values. In this regard, the problem with software-generated microarray gene signal intensities becomes more evident at the time of their fold change determination. This challenge is not only observed with the calculated fold changes, but also with the corresponding p-values, which may not be statistically significant (e.g., >0.05) and hence, relevant genes may be filtered out. However, with the normalization method described above, along with the utilization of reliable endogenous control genes, this can be corrected.

One important step taken before the analysis of gene expression is to restrict the search to a specific gene list or lists pertaining to a more focused field of interest (e.g., a particular disease-related list of genes). This helps in the manageability of the data set by restricting it to a much lower number of genes. The next step will be to filter the data according to the filtering parameters depicted in the previous section. That is, filtering by flags, leaving out those having compromised signals, and then filtering by expression, leaving out genes whose processed (background-subtracted) raw signal values are less than twice the background (S/N <2). Again, this is achieved by setting the processed signal’s lower limit to the equivalent of the technology’s mean background signal value. Once this selected group of genes is filtered, the next step is to manually scale the gene’s processed raw signals to the 75^th percentile, as described above. The resulting target gene values are then normalized by simple division using the combined value (i.e., mean) of the endogenous control or reference genes selected earlier, as follows: target gene value/endogenous control mean value, for each target gene. In this way, the values obtained can be used to compare control and treatment groups through fold change calculations (e.g., treated vs. control group) and the calculated p-values used to evaluate their statistical significance.

Conclusions

The ZDF rat is a proven model for the study of different comorbidities associated with type 2 diabetes. The results obtained in the present study demonstrate how use of a simple combination of software-generated and manual normalization methods can correct for potentially false results and aid in the selection of suitable endogenous control genes to be used in further gene expression determinations; in the present case, in the study of ZDF rat whole-blood samples. The expression of these genes showed no statistically significant differences between homozygous ZDF rats exhibiting clinically-relevant type 2 symptomatology and the heterozygous healthy lean controls, a characteristic which rendered them suitable. Importantly, the endogenous control genes that were found constitute a reliable platform for use in gene expression studies aiming to evaluate potentially novel therapeutic interventions for treatment of comorbidities and their progression in human populations with type 2 diabetes.

This method is especially useful when aimed to aid the software in cases of borderline results, where the expression and/or the fold change values are just beyond the pre-established set of acceptable parameters. In this regard, the difference between a gene with a p-value of 0.049 and one with a p-value of 0.051 is meaningless per se, as their true relevance is their biological significance. Hence, the use of endogenous control genes for the normalization of target genes assists in accomplishing the identification of potential biological significance in gene expression patterns. Moreover, this method should be applied every time in every array studied since, as noted earlier, differences in hybridization efficiencies along with changes in gene expression with different treatments and tissue-specific differential gene expression patterns are common occurrences.