caCORRECT2: Improving the accuracy and reliability of microarray data in the presence of artifacts

Artifact removal and replacement

Many artifacts were easily visible using the heat map function provided by caCORRECT, and one such example of mixed artifacts observed in published data is shown in Figure 1. The input and output of caCORRECT's artifact identification procedure are visualized, demonstrating the process by which continuous scoring of the quality of individual probes (dark color indicating poor quality in panel A) is combined to determine the binary classification of "artifact or not" (artifacts shown in bright red in panel B). Note that not every probe with a poor score was flagged as an artifact, as local areas of high variance are expected as a result of actual variance in gene expression among samples. For comparison, the artifact labels provided by Harshlighting are shown for both diffuse and compact artifacts in panel C and D respectively, with black and white indicating artifacts. This example supports the general trend that caCORRECT's automated artifact identification is conservative. It is our experience that human observers often consider a larger portion of a chip more objectionable than caCORRECT does. Harshlighting on the other hand, tends to be more aggressive in flagging regions of chips as artifactual. We chose conservative artifact identification over a more aggressive scheme because caCORRECT was designed to identify spatial artifacts upstream of the more specific model-based outlier detection employed by most probe summarization methods, such as RMA, PLIER, and MAS5.0. Generally speaking, probe summarization is the fitting of a regression model to the probe intensity data from a chip or set of chips. While individual methods differ on how the model is set up (notably how they handle information from mismatch probes), they are the same in that the model-fitting procedures allow for different weights or affinities for each probe's contribution to gene expression. During this model-fitting process, single outlier probes are easily detected and ignored because they do not fit the model well. caCORRECT is not designed to detect or remove such isolated, single probe defects, because RMA, MAS5.0 and PLIER are already very good at this task. With large artifacts that affect many probes, however, model-based outlier removal becomes much more difficult. To help users interpret the success or failure of the probe summarization (i.e. model-fitting) process, client software such as RMAExpress, and now caCORRECT, visualizes the residuals (difference, or error) of every probe intensity from the model collectively as a heat map. Close interpretation of these residuals reveals how well a particular model fits the data.

Figure 2 represents a case study using residual images to compare how caCORRECT and RMA react to the presence of a scratch on an Affymetrix Hu-6800 chip from the West dataset. We chose a data set based on this older array platform, because, on the Hu-6800 chip, all probes which contribute to a single gene expression value are actually arranged in contiguous regions on the microarray (as opposed to newer chip layouts which distribute related probes randomly). This property allows for easier visual interpretation of residual images. Panel A of Figure 2 shows the residuals produced by the caCORRECT regressive gene expression model before any attempt to identify artifacts. The blue (negative residual) regions that surround the main red scratch (positive residual) demonstrate the ambiguity that arises when more than one outlier is present in a probe set. Without prior knowledge, the regressive model cannot determine whether the artifact data are too high, or the non-artifact data are too low. caCORRECT, however, identifies the actual scratched region as artifact, and the surrounding data as non-artifact. With this knowledge, the regressive model fits only to the non-artifact data. The resulting residuals that are calculated after caCORRECT are shown in panel B of Figure 2. The white region surrounding the scratch (and redness of the scratch itself) is evidence that caCORRECT has successfully disambiguated the decision. The data within the scratch are to be ignored, and the resulting model should use only the surrounding probe data to estimate gene expression. Panel C of Figure 2 shows residuals to the model produced by RMAExpress alone. RMA includes its own outlier detection, but RMA's detection is not informed by spatial location, and thus cannot recognize a scratch as such. The output from RMA suggests that RMA does a fair job of disambiguating the scratch from the surrounding data, but some areas of blue still exist near the scratch, suggesting that the gene expression values for these probe sets have been overestimated. Finally, the residual image produced after the median data replacement of Harshlighting, shown in panel D, suggests that most if not all of some probe set calculations are based solely on the median-replaced data, with non-artifact data in surrounding regions being regarded as outlier data. These problems have been almost completely avoided by caCORRECT, however.

Effect of artifact aware quantile normalization on synthetic artifact data

Figure 3 illustrates the difference between standard quantile normalization and artifact-aware normalization. In this experiment the Schuetz et al. dataset has been modified with two large artifacts on one chip for illustrative purposes. As can be seen in the raw intensities (top panel), the synthetic artifacts (leftmost and rightmost modes in the red curve) cause a differently shaped distribution than that of a high quality chip (blue curve). Once standard quantile normalization is performed, a 'warping' of the high quality chip can be seen in the way that each of the distributions is now identically and incorrectly trimodal (middle panel). After caCORRECT artifact-aware normalization has been performed, the high quality chip returns to its natural distribution (bottom panel, blue curve). The chip with the artifacts still has two clear modes for its artifacts, but now the remaining non-artifactual data on the chip with the artifacts (central mode, red curve) has been properly normalized and aligned to the data from the high-quality chip (blue curve).

Effect of applied artifacts and preprocessing on gene expression

We monitored the effect of applied quality insults on gene expression using the two popular probe summarization methods of RMA and MAS5.0. Figure 4 demonstrates the correspondence of gene expression, before and after introduction of artifacts by an independent third party. Two phenomena are readily observable:

First, for the "black hole" artifacts which lower probe intensities on the microarray, the MAS5.0 algorithm has the tendency to call many of the genes 'absent', and report the gene expression abnormally low (Figure 4 panel C). caCORRECT is able to almost completely reverse this trend, and is able to help MAS5.0 produce appropriate gene expression values for most of these probe sets.

Second, for the "hot spot" artifacts that raise probe intensities on the microarray, the RMA algorithm has the tendency to underestimate gene expression and lose accuracy for the genes most highly expressed in the sample (Figure 4 panel B, red). This is presumably a result of the warping issues related to quantile normalization discussed in the previous section. This phenomenon also happens to low-expressing genes in RMA to a lesser extent for the black hole artifacts (Figure 4 panel D, red). Chips processed first with caCORRECT and then with RMA do not exhibit either of these warping behaviors (Figure 4, blue).

We then created our own synthetic insults in order to compare caCORRECT with Harshlighting for the ability to moderate effects of single artifacts on gene expression, as measured by the probe summarization methods RMA, PLIER, and MAS5.0. As a control to these common methods, which include normalization and outlier detection inherently, we also measured expression with our simple x _b,p,j = θ_p,ja _b,p + ε _b,p,j regression method "TAXY" detailed in the methods section, which does not contain special considerations for outlier data. Three conditions were tested for each method of gene expression measurement: (1) using original data that was not preprocessed; (2) using data that was processed with caCORRECT; and (3) using data that was processed with Harshlighting. We calculated the error caused by a given artifact as the average absolute difference, in the log domain, between the original gene expression and the expression measured after introduction of artifact. This measure is similar in concept to average relative error in gene expression in the linear domain.

Figure 5 shows the effect of artifact magnitude on gene expression. When the magnitude of the applied artifact was increased, the magnitude of error then increased for data that was not preprocessed, as well as for data processed with Harshlighting. In contrast, the magnitude of observed error remained constant, or even decreased, for data processed with caCORRECT. We expect that this is due to the phenomenon that as artifacts become more severe, they also become more obvious to caCORRECT. Even though harsher artifacts introduced more overall noise to the system, caCORRECT became better at identifying and removing them.

For RMA, caCORRECT outperformed both Harshlighting and unprocessed data for artifacts stronger than 1.5 fold. For PLIER and MAS5, however, caCORRECT outperformed unprocessed data only for intensity reducing artifacts. This suggests that caCORRECT is not suitable for helping MAS5.0 or PLIER to identify subtle artifacts, but that in the case of such subtle artifacts, caCORRECT does not reduce performance.

Figure 6 shows the effect of artifact size on gene expression. With the magnitude of artifact fixed, an increasing artifact size increased error for each gene expression and preprocessing method tested. For RMA, caCORRECT outperformed both Harshlighting and unprocessed data for artifacts larger than 25 probes across. For PLIER and MAS5, however, caCORRECT outperformed unprocessed data only for intensity-reducing artifacts. Overall, results suggest that the size, magnitude, and sign of an artifact all affect the final error in gene expression differently, depending critically on the gene expression method being used.

Identification of differentially expressed genes

Having established that caCORRECT could improve the accuracy of gene expression derived from microarrays in the presence of spatial artifacts, we set out to determine if this improved accuracy of gene expression would translate to improved efficiency in identifying candidate genes to serve as biomarkers for disease. Briefly, we first identified differentially expressed genes between subtypes of renal cell carcinoma from the Schuetz et al. microarray data, then we performed a small PCR pilot study to verify these findings, and finally we tested a much larger cohort of genes and samples via a PCR core facility. While results could not be directly quantitatively compared between experiments, trends in gene expression, such as ratios among biomarkers were generally preserved. Figure 7 shows the trends in gene expression for the two most reliable biomarkers identified during this process: NNMT and PRKAB1.

During the first PCR pilot study, we found only anecdotal evidence suggesting that using caCORRECT on real data could improve the reliability of biomarker selection. We believed that this could probably be attributed to the relatively high quality of the original chips in the Schuetz et al. study, and thus decided to artificially reduce the quality of the data in order to amplify any affect that caCORRECT may have. Figure 8 shows examples of our synthetic artifacts (Panels C and D) applied to the Schuetz et al. dataset as well as the unaltered versions of those chips (Panel A and B) as visualized by the post-caCORRECT residual images. While the synthetic artifacts may appear more visually stunning than the artifacts "naturally" found in this dataset, they are comparable with those found on other microarrays, such as the one shown covering the left hand side of the chip in Figure 1.

Figure 9 shows Receiver Operating Characteristic (ROC) curve analysis of microarray prediction (test prediction) versus PCR validation, which was defined as the gold standard for this purpose. For all cases, microarray data were generally predictive of follow-up PCR status. Results show that (1) caCORRECT was able to moderately improve biomarker selection power for typical microarray data (area under the curve increase from 0.777 to 0.786); and (2) caCORRECT was able to improve biomarker selection power for microarrays influenced by serious artifacts (area under the curve increase from 0.723 to 0.751). Importantly, caCORRECT had no undesired effects on the predictive power of clean data. Thus, caCORRECT should be suitable for datasets of unknown quality; in this setting, reliability would be expected to improve if the arrays contain significant artifact, while it would be unlikely to degrade if array artifacts are insignificant.

Affycomp

Cope et al. have provided a benchmark [20] by which to assess gene expression from two spike in datasets, one on the HG-U133A and the other on the HG-U95A platform [8]. Their spike in arrays were all prepared using the same cRNA stock mixture, with the exception of 16 transcripts which were added in concentrations from 0 to 1024 picoMolar. These 16 transcripts were applied (in triplicate) in a cyclic Latin square design. As a result, the fold change between any two arrays for these 16 transcripts is known, while the fold change for all other transcripts is expected to be 1 (no change). The affycomp package, available from http://bioconductor.org/, provides many statistics and figures of merit with which to compare the accuracy of gene expression in the context of detecting known fold change between pairs of arrays. Affycomp provides a uniquely public way to assess caCORRECT without the need to synthetically worsen the data. Figure 10, panels A and B, show caCORRECT's performance on the HG-U95A data set as a plot of false positive rate versus true positive rate (ROC curve). RMA after caCORRECT was the best overall performer, as measured by Area Under the Curve (AUC) for less than 100 false positives (0.830 for RMA after caCORRECT, 0.821 for RMA alone, and 0.818 for RMA after Harshlighting). AUC gains from caCORRECT were most noticeable for the TAXY method, which does not provide any additional outlier detection like that inherent in RMA, MAS5.0 and PLIER (0.760 for TAXY after caCORRECT, 0.740 for TAXY alone, and 0.719 for TAXY after Harshlighting). Importantly, the poorest-quality chip (Figure 10, panel C) was thrown out by the original authors as part of the canonical affycomp benchmark-- a decision which was recently justified quantitatively [21]. Still, examples of chips on which caCORRECT has detected obvious artifacts (Figure 10, panel D) remain in the HG-U95 spike in data. For the HG-U133A spike in dataset, which had less detected artifacts, improvements due to caCORRECT were more modest

It is almost never a good idea to run caCORRECT on batches with 3 or fewer chips, although it is hard to imagine a reliable microarray experiment so small. Affycomp data, however, followed a rarely achieved design in which there are 3 technical replicates of each sample. Arteaga-Salas et al. were able to apply their method for correcting small batches of replicates by splitting up the affycomp data set [17]. They have previously reported performance as the fraction of spike in genes from affycomp HG-U133A to have their RMA-calculated fold-change rank improved or worsened as a result of correction. Processing batches of 3 technical replicates at a time, they report 45.70% improved ranks, while caCORRECT, processing the entire dataset at once, improved ranks of 10.81% of spiked in genes. However, they report 38.60% worsened ranks, while caCORRECT worsened the ranks of only 9.52% spiked in genes. This result is consistent with caCORRECT's conservative design, but suggests that Arteaga-Salas's correction method may also be appropriate for experimental designs with available technical replicates for every sample.

Variance scoring

The cornerstone of caCORRECT's outlier detection is the concept of variance scoring, which is a description of each probe's tendency to be an outlier. Calculation of this score, h, is similar to conducting a t-test for whether or not the observed probe intensity for a given chip belongs to the observed distribution of probe intensities for all other chips in the dataset. A key feature of caCORRECT is that this distribution is estimated and updated multiple times during the course of a single caCORRECT session. Because of this dynamic updating, it is possible to identify subtle artifacts or pardon false artifacts that may have been misdiagnosed initially. Please refer to File01_supplement.pdf, "Supplemental Methods", for a more detailed description of how the variance score is calculated.

Artifact segmentation

Once the variance statistic, h, has been calculated for each probe on each chip, false-color heat maps of h, showing probes in their original spatial layout, are generated to display regions of high noise. For a good quality microarray chip, h will represent biological variation in RNA expression for the sample. In this case, h will be distributed independently and nearly-uniform in magnitude throughout the chip. More commonly, however, protocols do not achieve uniform hybridization due to uneven drying, formation of salt streaks, scratching of the microarray surface due to contact with skin or dust, miscalculated hybridization times, or failure to control environmental variables such as ozone [24]. All of these most common mistakes result in localized regions of large h (artifacts) on the heat map.

caCORRECT uses a simple sliding window method to flag probes that meet two conditions: (1) they exist in regions of other high h-scoring probes, and (2) they have high h scores themselves. These two conditions ensure that most of the obvious artifacts are caught, but that most of the naturally occurring biological variance goes unnoticed. Because the intended platform for caCORRECT is a web-based grid service, artifact identification has been streamlined for speed and memory efficiency. More computationally intense methods such as active contours, PDE-based methods, or shape matching have been excluded in favor of a quick marching window algorithm that seems to work well for a wide range of data. To remove any global chip effects that arise from sample preparation or amplification, normalization is performed as described in the following section.

Artifact-aware normalization

Quantile normalization reduces noise in microarray experiment replicates by forcing the intensity distribution of each chip to be identical [Bolstad B: Probe level quantile normalization of high density oligonucleotide array data, 2001]. The critical assumption behind quantile normalization is that for large genome-wide studies such as microarrays, the number of genes that are invariant to the experimental variables far outnumbers the number of biomarkers-- genes that respond to or predict experimental variables. Quantile normalization is generally good for the microarray problem, where the distributions are poorly defined and parametric methods such as median centering or Z-score normalization have their underlying assumptions violated. The power of quantile normalization comes with a major caveat. If the probe intensities of the chips are not distributed similarly, the algorithm will indiscriminately warp all the distributions to be the same, including any that may have been correct initially. Fortunately, it is a reasonable assumption that high-quality microarray data from a single source on a single platform follow the same distribution. Unfortunately, this high quality assumption is not valid for much real-world data, where chip artifacts can significantly alter the distribution of intensities on a chip. One bad chip can warp the others when quantile normalization is performed, thus compromising the reproducibility of the entire data set. A way to alleviate this problem is to identify artifacts before quantile normalization, and set them aside temporarily. In theory, perfect knowledge of artifacts would allow for perfect correction. This process is called "artifact-aware quantile normalization." caCORRECT uses four iterations of artifact-aware normalization and artifact identification in order to achieve a near steady-state normalization result with a relatively small amount of computation time.

To illustrate the invasive effect that artifacts can have on a dataset when quantile normalization is performed, synthetic microarray data were generated in the following manner: Six high-quality chips from the Schuetz et al. dataset [25] were chosen, one of which was set aside to receive artifacts. One third of the selected chip was modified by a multiplicative factor of 0.5, representing a low-intensity artifact. A different third of the selected chip was modified by a multiplicative factor of 10, representing a high-intensity artifact. These six chips were then processed using caCORRECT, and the probe intensities were monitored for warping at each intermediate step.

Artifact replacement and probe intensity model

Once artifacts have been identified, a clean dataset is generated with the artifactual data appropriately replaced. The current version of caCORRECT uses a data imputation that is mathematically similar to the model used by RMA and others [5, 9, 10]. Notably, so-called perfect-match and mismatch probes are treated identically, i.e. we do not use PM-MM. In this scheme, the artifact-flagged probe level data is replaced with the best-fit estimate for that probe, given the model and data from non-artifactual probes on all of the chips being processed.

where x _b,p,j is the observed intensity for the b ^th probe in the p ^th probe set on the j ^th chip, θ _p,j is the gene expression term, a _b,p is the probe affinity term, and ε _b,p,j is the additive error term.

The solution which satisfies the above condition can be derived from the singular value decomposition (SVD) of X _p . Here, the SVD of X _p is given in the form of X _p = USV^T, such that,$U\in {\Re }^{N\times N}$, $S\in {\Re }^{N\times B_{\mathsf{\text{p}}}}$ and $V\in {\Re }^{B_{\mathsf{\text{p}}}\times B_{\mathsf{\text{p}}}}$. If the largest singular value in S, s₁, is arranged as the first diagonal element of S, then θ _p is s₁ times the first column of U and $a_p^{\mathsf{\text{T}}}$ is the first column of V.

We chose to introduce the additional constraint that the geometric mean of the lumped probe affinity terms a _b,p equals one. The number one is arbitrary here, but it allows the convenient interpretation that the values of gene expression, θ _p,j , are on the same scale as the probe intensities, x _b,p,j . To satisfy this constraint, the earlier solution for a _p and θ _p can simply be scaled by respective multiplication and division.

Imputation of artifact values

This procedure of replacing values in X _p with values from θ _p a _p has the effect of reducing corresponding values of ε _p to zero, and thus has the property of never increasing the Frobenius norm of ε _p . Since step 1 is a global minimization of the Frobenius norm of ε _p given X _p , and step 2 alters X _p in a way that can only further decrease this error, the entire procedure is guaranteed to converge due to the non-increasing nature of the error function, which is naturally bounded by zero. Troyanskaya et al. have previously used a similar "SVDImpute" procedure for missing gene expression data [22]. Please see additional file 1, "Supplemental Methods", for further details.

Datasets and synthetic artifact generation

In order to quantify the ability of caCORRECT to improve data quality, we altered public microarray datasets with a variety of randomized synthetic artifacts, and then processed the altered data with caCORRECT or Harshlighting. Datasets were generated from a variety of clinical cancer studies using different microarray platforms. To date, the caCORRECT website has been tested with data from 18 different Affymetrix platforms, but the results of our synthetic artifact analysis which are presented here are limited to our in-depth study of two key data sets. Three separate experiments were performed involving application of synthetic artifacts.

Third party artifacts on Hess et al. data

First, we obtained a large data set that was originally generated by Hess et al. [26] in the study of breast cancer, and then later used as part of MAQC phase-II study on classifier performance. This original data set consisted of 130 high-quality samples assayed on the Affymetrix HG-U133A platform. The Hess et al. study divided samples into training (n = 81) and validation (n = 49) sets and we retain this distinction. Chips in the validation set were selected for synthetic artifact manipulation by an independent team lead by Wendell Jones of Expression Analysis, (previously unaffiliated with caCORRECT). Two types of artifacts were investigated here: (1) a "black hole" artifact in which an elliptical region of the microarray had probe intensities lowered severely, and (2) a "hot spot" artifact in which an elliptical region of the microarray had probe intensities raised severely. Twenty digitally altered copies of each of the original 49 chips were prepared as follows: A single artifact with random orientation and location was applied to each chip (ten chips received "black holes", and ten received "hot spots"). For each of the altered chips, gene expressions were calculated both before and after caCORRECT's complete artifact detection and value imputation. Expression data for all probes were determined both using MAS5 in Expression Console and the R implementation of RMA. Each of the estimated gene expressions from the altered chips was then compared to the "true" gene expression values obtained from the respective original, unaltered chip to yield an error value representing the deleterious effect of the artifact on gene expression estimation. The errors for each probe set (22283), each chip (49), and each artifact replicate of a given type (10) were then pooled together to form distributions of errors of size n = 10,918,670. Eight such distributions were created in total, representing the combinations of two gene expression methods, two artifact types and either unprocessed data, or for data cleaned with caCORRECT.

Artifacts generated on data from Hess et al

A common criticism of our synthetic artifact work is that the size and severity of tested artifacts (for example those provided by Jones' team) are rarely observed in practice. While we have encountered many chips plagued by large, severe artifacts, (see http://arraywiki.bme.gatech.edu/index.php/Hall_of_Fame for some examples [27]), we set out here to investigate how the magnitude and size of artifacts may affect our previous conclusions as to the usefulness of caCORRECT. Thus, we have created a second set of synthetic artifact data specifically for this purpose by our own application of artifacts to the same breast cancer dataset from Hess et al that Jones' team used. Only the 15 highest quality arrays (via visual inspection of heat map images) from this dataset were used in order to both speed up computation as well as to more precisely isolate the effects of our synthetic artifacts. For a variety of sizes and magnitudes, circular regions of the array were altered multiplicatively. Care was taken so that no more than 1/2 of the radius of the circular insult appeared off of the chip. The final gene expression obtained from the altered array was then compared to the gene expression obtained from the original unaltered array, and the average of the relative error for each probe set on the array was stored. For each pair of size and magnitude, this process was repeated a total of 3 times.

Artifacts observed and generated on data from Schuetz et al

Finally, a realistic mixture of less-severe artifacts were applied to the Schuetz et al. [25] dataset in order to monitor the effect of typical artifacts on differential gene finding, and the ability of caCORRECT to ameliorate these effects. This data set consisted of 20 Renal Cell Carcinoma (RCC) samples assayed on the Affymetrix HG-Focus Platform by Schuetz et al. Samples were classified by tumor subtype: Clear Cell (CC), Oncocytoma (ONC), and Chromophobe (CHR). For biomarker selection purposes, samples were combined into two classes: seven CHR or ONC versus thirteen CC.

Artifacts observed on data from West et al

The dataset which was used to showcase real-world artifact removal and replacement was a set of 49 Hu-6800 Affymetrix microarrays from the study by West et al. This study investigated Estrogen Receptor and lymph node metastatic status [28]. This data set is chosen because it used an older version of Affymetrix chip in which the properties of artifacts are more easily visualized than on modern chips.

PCR Validation

To determine the effect that caCORRECT had on the ability to correctly identify biomarkers of disease from microarray data, a panel of 96 genes of interest for RCC was assembled for PCR study in two phases. These genes were identified from a combination of genes previously identified in the literature as well as a set of genes whose biomarker status was disagreed upon between the caCORRECT and non-caCORRECT versions of the Schuetz et al. data sets using omniBiomarker (http://omniBiomarker.bme.gatech.edu). All PCR analysis was performed on independent patient tissue samples with respect to those used for the microarray analysis.

Gene expression was assessed by quantitative RT-PCR, using total RNA from fixed tissues of 17 clear cell and 7 chromophobe RCC patients. Duplicate experiments were performed according to published protocols with minor modifications: Histological sections were deparaffinized with ethanol and xylene, and cells of interest were microdissected with a sterile scalpel. Tissues were digested in buffer containing proteinase K at 60°C overnight. RNA was extracted with phenol/chloroform, and genomic DNA was removed with DNase. RNA quality and quantity were assessed with a Bioanalyzer (Agilent Technologies). Up to 3 μg of RNA was used for first strand cDNA synthesis with Superscript III (Invitrogen). PCR was performed with a custom-designed Taqman Low Density Array (LDA, Applied Biosystems) in a 96-well microfluidic card format, using the ABI PRISM 7900HT Sequence Detection System (high-throughput real-time PCR system). Gene expression data were normalized relative to the geometric mean of two housekeeping genes (18S, ACTB). LDA runs were analyzed by using Relative Quantification (RQ) Manager (Applied Biosystems) software.

Relative normalized PCR gene expression was compared in renal tumor subtypes. Genes were declared as being "validated by PCR" if they had an average fold change between classes of magnitude greater than 2, corresponding to an average C_t (threshold cycle) value difference of 1.