Using control genes to correct for unwanted variation in microarray data

Abstract

Microarray expression studies suffer from the problem of batch effects and other unwanted variation. Many methods have been proposed to adjust microarray data to mitigate the problems of unwanted variation. Several of these methods rely on factor analysis to infer the unwanted variation from the data. A central problem with this approach is the difficulty in discerning the unwanted variation from the biological variation that is of interest to the researcher. We present a new method, intended for use in differential expression studies, that attempts to overcome this problem by restricting the factor analysis to negative control genes. Negative control genes are genes known a priori not to be differentially expressed with respect to the biological factor of interest. Variation in the expression levels of these genes can therefore be assumed to be unwanted variation. We name this method "Remove Unwanted Variation, 2-step" (RUV-2). We discuss various techniques for assessing the performance of an adjustment method and compare the performance of RUV-2 with that of other commonly used adjustment methods such as Combat and Surrogate Variable Analysis (SVA). We present several example studies, each concerning genes differentially expressed with respect to gender in the brain and find that RUV-2 performs as well or better than other methods. Finally, we discuss the possibility of adapting RUV-2 for use in studies not concerned with differential expression and conclude that there may be promise but substantial challenges remain.

Fig 1.

Gender study RLE plots at different stages of preprocessing. From left to right: No preprocessing; BG/QN done separately for each platform type; BG/QN followed by a final LS across all chips. Coloring: red—site A, HG-U95A; yellow—site A, HG-U95Av2; black—site B, HG-U95A; gray—site B, HG-U95Av2; and cyan—site C, HG-U95Av2. Note: The scale on the y-axis is different for these RLE plots than for all other RLE plots in this paper.

Fig 2.

Gender study p-value histograms and RLE plots before and after adjustment. Histogram breakpoints are at 0.001, 0.01, 0.05 and 0.1, 0.2, 0.3, etc. Data were fully preprocessed (BG + QN + LS). The factors were computed by SVD on the HK genes. P-values were computed using Limma.

Fig 3.

Comparison of performance of different factor analysis methods in the gender study. The number of X/Y genes discovered is plotted as a function of k. Genes were ranked by p-value; results are shown for the number of X/Y genes ranked in the top 20 (plus), top 40 (triangle), and top 60 (circle). $k=0$ corresponds to no adjustment. Data were preprocessed (BG + NM + LS). Principal components were computed using the HK genes. P-values were computed using Limma.

Figure 4

Comparison of results for HK genes and Affymentrix spike-in controls in the gender study. For the RLE plots and p-value histograms, $k=10$. Factors were computed by SVD. P-values were computed using Limma. Note that there are only 33 spike-in controls, so adjustments with $k>33$ are undefined for the spike-in case. We truncate