Preferred analysis methods for Affymetrix GeneChips revealed by a wholly defined control dataset

Abstract

Background: As more methods are developed to analyze RNA-profiling data, assessing their performance using control datasets becomes increasingly important.

Results: We present a 'spike-in' experiment for Affymetrix GeneChips that provides a defined dataset of 3,860 RNA species, which we use to evaluate analysis options for identifying differentially expressed genes. The experimental design incorporates two novel features. First, to obtain accurate estimates of false-positive and false-negative rates, 100-200 RNAs are spiked in at each fold-change level of interest, ranging from 1.2 to 4-fold. Second, instead of using an uncharacterized background RNA sample, a set of 2,551 RNA species is used as the constant (1x) set, allowing us to know whether any given probe set is truly present or absent. Application of a large number of analysis methods to this dataset reveals clear variation in their ability to identify differentially expressed genes. False-negative and false-positive rates are minimized when the following options are chosen: subtracting nonspecific signal from the PM probe intensities; performing an intensity-dependent normalization at the probe set level; and incorporating a signal intensity-dependent standard deviation in the test statistic.

Conclusions: A best-route combination of analysis methods is presented that allows detection of approximately 70% of true positives before reaching a 10% false-discovery rate. We highlight areas in need of improvement, including better estimate of false-discovery rates and decreased false-negative rates.

Figure 1

Schematic depiction of the experimental protocol.

Figure 2

Signal of individual probes and dependence on present versus absent RNA molecules. (a, b) Plot of probe-pair signals for the three C chips, highlighting (a) the empty probe pairs or (b) the present probe pairs in green. (c) Receiver-operator characteristic (ROC) curves at the probe-pair level for several absent/present metrics. The metric (PM - MM)/(PM + MM) gives the same result as the green curve. (d)Receiver-operator characteristic curves at the probe-set level for several absent/present metrics combined using the Wilcoxon rank sum test.

Figure 3

The set of options that were investigated using Bioconductor's affy package. The choices that optimize the detection of DEGs are circled in red. Broken circles indicate choices that are slightly suboptimal but still rank within the top 10 datasets.

Figure 4

The dependence of log fold change on signal intensity (M versus A plots). (a)M versus A plot before the second normalization step and (b) after a loess fit at the probe set level. FC in the key denotes the spiked-in fold change value.

Figure 5

Correlation of observed with actual fold changes for a representative expression summary dataset (Additional data file 2, using dataset 9e.b). (a) The fold change for each probe set with spiked-in target RNA is depicted as a cross. Empty probe sets are not shown. For each actual fold-change level (on the x axis), a boxplot shows the distribution of the corresponding observed fold changes. A linear fit of the data is shown in cyan. Fit parameters: R² = 0.508; slope = 0.505; y-intercept = -0.061. (b-d) Increasingly more of the low-intensity probe sets are filtered out of the plot. All probe sets are ranked according to average signal level, and those in the lowest 25th (b), 50th (c), or 75th (d) percentile of signal level are eliminated from (a). Fit parameters: (b) R² = 0.870; slope = 0.546; y-intercept = -0.008; (c) R² = 0.895; slope = 0.517; y-intercept = -0.015; (d) R² = 0.906; slope = 0.457; y-intercept = -0.017.

Figure 6

Comparison of three t-statistic variants. (a)ROC curves for a particular expression summary dataset, using the different t-statistics. Location of false positives and false negatives are shown for the (b) CyberT, (c) SAM, and (d) basic t-statistic when considering the top 1,000 probe sets as positive DEG calls.

Figure 7

ROC curves for all expression summary datasets. The curves are color-coded to highlight how the ability to detect differential expression is dependent on the different options at each step of analysis, using the CyberT regularized t-statistic metric. (a) All 152 expression summary datasets are represented here, with the different colors depicting whether the second loess normalization step at the probe set level was performed. In general, the second loess normalization (blue) improves the detection of true DEGs. (b-f)To decrease clutter, only the 76 expression summary datasets involving the second normalization step are shown. (b) When comparing the two background correction methods, the MAS algorithm is superior to the RMA algorithm. (c) The various probe-level normalization methods do not show great differences between each other. (d) Among the different PM-correction options, using the method in MAS 5.0 clearly is the most successful. (e) Various robust estimators were examined, revealing that the median polish method is the most sensitive (with MAS 5.0's Tukey Biweight a close second). (f) Depiction (in blue and orange) of the 10 datasets which maximize detection of truly differentially expressed genes, while minimizing false positives. These datasets are generated using the options circled in Figure 3. MAS 5.0, with the inclusion of the second loess normalization step, falls within these top 10.

Figure 8

The accuracy of false discovery rate estimates (q-values). The top 10 expression summary datasets (named 9a-9e, 10a-10e in Additional data file 2) were combined to generate a composite statistic, which was used to rank genes based on the robustness of their significance over the 10 datasets. (a) The composite statistic performs as well as the best summary dataset in terms of sensitivity and specificity. (b) In addition, permutation tests carried out using this composite statistic yield q-value estimates which are more accurate than any of the 10 component datasets, although still lower than the true false-discovery rate.

Figure 9

DEG detection sensitivity and specificity as a function of spiked-in fold change level. (a, b)ROC curves using the composite statistic, and different definitions of the true-positive probe sets (criteria given in the legends; FC, spiked-in fold change). The true negatives remain the same for all curves (the probe sets which were not spiked in, or were spiked in at 1x).