Impact of gene expression data pre-processing on expression quantitative trait locus mapping

Evidence for heritability of mRNA levels has been observed in several organisms like the mouse, yeast, and human [1–3]. Therefore, it is possible to consider transcript levels measured using DNA microarrays as quantitative traits and localize the genes controlling them by quantitative trait locus (QTL) analyses. However, one of the challenges is to remove optical noise and take into account nonspecific hybridization in the microarray experiments. In addition, in the Affymetrix system, each gene is represented by 11–20 Perfect Match (PM) and Mismatch (MM) pairs of probes, each probing a different region of the mRNA transcript, typically within 600 base pairs of the 3' end [4]. It is also important to find a way to combine the 11–20 probe pair intensities for a given gene into a single measure. Several data pre-processing methods are now available. Of several such proposed methods, three have commonly been employed: the Robust Multiarray Average (RMA), GeneChip RMA (GCRMA), and the Microarray Analysis Suite 5 (MAS5) methods. A detailed description of these methods can be found elsewhere [5–7]. These three methods convert probe-level data to expression values through the following sequence: 1) background correction, 2) normalization, and 3) summarizing the probe set values into one expression measure. Regarding the background correction step, MAS5 uses the MM probes to adjust the PM probes for probe-specific nonspecific binding [5]. On the other hand, RMA ignores the MM intensities and performs a global background correction [6]. GCRMA employs a hybrid approach and uses the probe sequence information released by Affymetrix to compute an affinity measure and describe background noise. For Step 2, RMA and GCRMA use the same normalization method based upon intensity quantiles, which imposes the same empirical distribution of intensity to each array. In MAS5, a baseline array is chosen and all the other arrays are scaled to have the same mean intensity as this array. Finally, RMA and GCRMA use the same probe summary method, based on a robust linear model. In MAS5, the probe summary is based on a robust average method.

Different DNA chip pre-processing methods have been shown to influence measures of gene expression [6, 8, 9]. Therefore, the choice of a pre-processing method significantly affects linkage results obtained in genomic analyses [10]. However, whether or not only concordant results should be further investigated is still debated (see Petretto et al. [11] and Chesler et al. [12] in their reply to Williams et al. [10]). In this paper, we focus on comparing the degree of concordance between pre-processing methods, according to the characteristics of the traits selected for linkage analysis.

The study is based on all 14 three-generation CEPH (Centre d'Etude du Polymorphisme Humain) Utah families. Gene expression levels in lymphoblastoid cells of 194 individuals have been obtained using the Affymetrix Human Focus Arrays that contain probes for 8792 transcripts. Details regarding the microarray experiments are given in Morley et al. [3].

As previously suggested by other studies [6, 8, 10], our results confirm that pre-processing methods may also affect linkage outcomes. However, this impact depends on the way traits were selected for genetic analysis. Choosing the traits on the basis of a high heritability value led to a minimum discrepancy between methods. Conversely, discrepancies were more important in the group of non-expressed traits or in the group of variable traits chosen at random (note that in this last group, heritability ranges from 0% to 50%). Furthermore, we noticed that the three pre-processing methods agree much better for cis-acting than for trans-acting regulators.

Several factors may explain partly why the three methods produce different results. First, as already stated in the Background, the underlying models converting probe level data to expression values are different from one method to another. Although the normalization step is not the same for the three methods, previous work has shown that these differences have little effect relative to that of the background correction, which entails a variance/bias trade-off. Especially, it has been shown that background correction decreases the bias but that naïve background correction procedures, such as MAS5 and RMA, increase the variance [8]. GCRMA is supposed to provide a good balance between accuracy and precision by doing adequate non-specific binding correction. Our study suggests that the large impact of the background correction also applies to eQTL mapping results. Indeed, RMA and GCRMA differ only by the background correction step, yet their concordance rates were not particularly high. Interestingly, Irizarry et al. [8] observed that differences in precision between RMA and GCRMA were higher in the case of genes with low expression, with GCRMA giving the smallest bias. We also observed a greater discordance rate between RMA and GCRMA for the set of non-expressed genes.

Another potential factor for explaining differences between methods is departure from normality of the phenotypic distribution, especially when using a variance-component approach. We found that, in general, GCRMA led to the highest rate of traits failing the Shapiro normality test. Among expressed genes, these rates were 5.4, 0.8, and 1.7% for GCRMA, RMA, and MAS5, respectively (using a Bonferroni correction for multiple testing at level 5%). These rates might explain the large number of linkage signals observed for a few traits with GCRMA. However, it seems unlikely that conflicting eQTL mapping results are mainly due to differences in the gene expression distributions per se.

In conclusion, the true genetic determinants of the studied traits in the GAW Problem 1 data are unknown, preventing us from drawing definite conclusions on the best and more robust pre-processing method. Further, in the context of a genome-scan, high agreement rates across experiments are not expected because most of the linkage signals are likely to be false positives. It is unclear whether it would be sound to use several pre-processing methods in a systematic manner. Such guidelines were proposed recently but remain controversial [10–12]. In our study, we found very poor agreement between pre-processing methods in the set of the non-expressed but most variable genes (i.e., our Set 4), suggesting that filtering genes on their detectable presence in the tissue of interest is also an important step. To filter genes not only on their variability but also on their presence in the tissue analyzed is one of our main messages.