Removing batch effects for prediction problems with frozen surrogate variable analysis

Batch effects are responsible for the failure of promising genomic prognos- tic signatures, major ambiguities in published genomic results, and retractions of widely-publicized findings. Batch effect corrections have been developed to re- move these artifacts, but they are designed to be used in population studies. But genomic technologies are beginning to be used in clinical applications where sam- ples are analyzed one at a time for diagnostic, prognostic, and predictive applica- tions. There are currently no batch correction methods that have been developed specifically for prediction. In this paper, we propose an new method called frozen surrogate variable analysis (fSVA) that borrows strength from a training set for individual sample batch correction. We show that fSVA improves prediction ac- curacy in simulations and in public genomic studies. fSVA is available as part of the sva Bioconductor package.

surements made by microarrays can be used to classify patients into predictive, prognostic, or diagnostic groups. Despite the incredible clinical promise of these technologies there have only been a few signatures that have successfully been translated into the clinic.
One of the reasons for the relatively low rate of success is the impact of unmeasured technological or biological confounders. These artifacts are collectively referred to as "batch effects" because the processing date, or batch, is the most commonly measured surrogate for these unmeasured variables in genomic studies [23,12,28]. The umbrella term batch effects also refers to any unmeasured variables that can vary from experiment to experiment, ranging from the technician who performs the experiment to the temperature and ozone levels that day [14,9].
Batch effects are responsible for the failure of promising genomic prognostic signatures [2,3], major ambiguities in published genomic results [25,1], and retractions of widely-publicized findings [24,13]. In many experiments, the signal from these unmeasured confounders is larger than the biological signal of interest [18]. But the impact of batch effects on prediction problems has only recently been demonstrated [22,20]. Batch effects were also recognized as a significant hurdle in the development of personalized genomic biomarkers in the Institute of Medicine's report on clinical genomics [21].
While a number of methods have been developed for removing batch effects in population-based genomic studies [12,11,16,18,28], there is currently no method for removing batch effects for prediction problems. There are two key differences between population level corrections and corrections designed for prediction problems. First, population level corrections assume that the biological groups of interest are known in advance. In prediction problems, the goal is to predict the biological group. Second, in prediction problems, new samples are observed one at a time, so the surrogate batch variable will have a unique and unknown value for each sample.
Here we propose frozen Surrogate Variable Analysis (fSVA) as a method for batch correction in prediction problems. fSVA borrows strength from a reference database to address the challenges unique to batch correction for prediction. The fSVA approach has two main components. First, surrogate variable analysis (SVA) is used to correct for batch effects in the training database. Any standard classification algorithm can then be applied to build a classifier based on this clean training data set. Second, probability weights and coefficients estimated on the training database are used to remove batch effects in new samples. The classifier trained on the clean database can then be applied to these cleaned samples for prediction.
We show with simulated data that the fSVA approach leads to substantial improvement in predictive accuracy when unmeasured variables are correlated with biological outcomes. We also apply fSVA to multiple publicly available microarray data sets and show improvements in prediction accuracy after correcting for batch. The methods developed in this paper have been implemented in the freely available sva Bioconductor package.
2 Frozen Surrogate Variable Analysis methodology 2.1 Removing batch effects from the training set The first step in batch correction for prediction problems is to remove batch effects from the training set. In the training set, the biological groups are known. This setting is similar to the population genomics setting and we can use a model for gene expression data originally developed for population correction of unmeasured confounders. If there are m measured features and n samples in the training set, we let X m×n be the matrix of feature data, where x i j is the value of feature i for sample j. For convenience, we will refer to X as the expression matrix for the remainder of the paper. However, our methods can be generally applied to any set of features, including measures of protein abundance, gene expression, or DNA methylation.
We propose a linear model for the relationship between the expression levels and the outcome of interest y j : The s k (·) are a set of basis functions parameterizing the relationship between expression and outcome. If the prediction problem is two-class, then p 1 = 1 and s 1 (y j ) = 1(y j = 1) is an indicator function that sample j belongs to class one. In a multi-class prediction problem k > 1 and the s k (·) may represent a factor model for each class. In matrix form this model can be written [16,17].
where S p 1 ×n is a model matrix of p 1 biological variables of interest for the n samples, B m×p 1 are the coefficients for these variables, and E m×n is the matrix of errors.
In genomic studies, the error term E is not independent across samples [12,16,17,28,10,18,11]. That is, there is still correlation between rows of E after accounting for the model S. The correlation is due to unmeasured and unwanted factors such as batch. We can modify model 1 to account for these measured biological factors and unmeasured biological and non-biological factors: where G p 2 ×n is a p 2 × n random matrix, called a dependence kernel [17] that parameterizes the effect of unmeasured confounders, Γ m×p 2 is the m × p 2 matrix of coefficients for G, and U m×n is the m × n matrix of independent measurement errors. We previously demonstrated that such a decomposition of the variance exists under general conditions typically satisfied in population genomic experiments [17].
In the training set, the biological classes are known, so S is known and fixed. But the matrices B, Γ and G must be estimated. fSVA first performs surrogate variable analysis (SVA) on the training database in order to identify surrogates for batch effects in the training samples. The training set can be "cleaned" of batch effects by regressing the effect of the surrogate variables out of the data for each feature. Any classification algorithm can then be developed on the basis of the clean training data set.
SVA is an iterative algorithm that alternates between two steps. First SVA estimates the probabilities π iγ = Pr(γ i· = 0|X, S,Ĝ), π ib = Pr(b i· = 0|γ i· = 0, X, S,Ĝ) using an empirical Bayes' estimation procedure [17,7,26]. These probabilities are then combined to define an estimate of the probability that a gene is associated with unmeasured confounders, but not with the group outcome The second step of the SVA algorithm weighs each row of the expression matrix X by the corresponding probability weightpi iw and performs a singular value decomposition of the weighted matrix. Lettingŵ ii =π iw the decomposition can be writtenŴ X = UDV T . After iterating between these two steps, the first p 2 weighted left singular vectors of X are used as estimates of G. An estimate of p 2 can be obtained either through permutation [4] or asymptotic [15] approaches.
Once estimatesĜ have been obtained, it is possible to fit the regression model in equation 2 using standard least squares. The result are estimates for the coefficientsB andΓ. Batch effects can be removed from the training set by settingX clean = X −ΓĜ.
Any standard prediction algorithm can then be applied toX clean to develop a classifier based on batch-free genomic data. The result is a prediction function f (X clean · j ) that predicts the outcome variable y j based on the clean expression matrix. 3 Fast fSVA methodology fSVA requires that a new singular value decomposition be applied to the augmented expression matrix once for each new sample. Although this is somewhat computationally intensive, in typical personalized medicine applications, sample collection and processing will be spread over a long period of time. In this setting, computational time is not of critical concern. However, for evaluating the fSVA methodology or developing new classifiers using cross-validation, it is important to be able to quickly calculate clean expression values for test samples.

Removing batch effects from new samples
We propose an approximate fSVA algorithm that greatly reduces computing time by performing a streaming singular value decomposition [29,30]. The basic idea be- Estimates obtained using this approximate algorithm are not identical to those obtained using the exact fSVA algorithm. The projection matrix used in the approximation, P T , is calculated using only the samples in the training set. However, there is only a one-sample difference between the projection calculated in the training set and the projection that would be obtained with exact fSVA. As the training set size grows, the approximation is closer and closer to the answer that would be obtained from the exact algorithm. For smaller databases, there is less computational burden in calculating the exact estimates. However, for large training sets, the computational savings can be dramatic, as described in the simulation below.
We performed a simulation examining the benefit of fSVA in prediction problems. In order to do this, we simulated data using equation 2 under different distributions of each parameter. We also discretized the probability weights π iγ and π ib , and varied the distribution of these probability weights (Table 4). We crafted these simulations to mimic scenarios with a subtle outcome and a strong batch effect, which is frequently the case in genomic data.

Parameter Distributions Affected Features
40% affected by batch and outcome 40% affected by batch and outcome 50% affected by batch and outcome Table 1: Specifications for the three simulation scenarios used to show the performance of fSVA. We performed three simulations under slightly different parameterizations to show the effectiveness of fSVA in improving prediction accuracy. Parameters from equation 2 were simulated using the distributions specified in this table. Additionally, the percentage of features in the simulation affected by batch, outcome, or both are as indicated in this table. Results from these simulations can be found in Figure 1.
We also specified that each data set have two batches and two outcomes, each equally distributed. We varied the amount of confounding between batch and outcome in the simulated database, from a Pearson's correlation of 0 to a correlation of over 0.90.
For the simulated new samples, the batch and outcome were always uncorrelated.
We simulated 100 database samples and 100 new samples using the parameters described above. Each sample had 10000 features. As a control, for each iteration in addition to performing fSVA correction, we performed SVA correction on the simulated database alone, and also performed prediction with no batch correction on the simulated database or new samples.
To quantify the effect that fSVA had on prediction, we performed exact fSVA as described above on the simulated database and new samples. We then performed Prediction Analysis of Microarrays (PAM), a commonly used method for classifying microarrays [27]. The PAM prediction model was built on the SVA-corrected database, and then used to predict the outcomes on the fSVA-corrected new samples. Each simulation was repeated 100 times for robustness.
We found that in general the prediction accuracy measures for different iterations of a simulation varied highly, but the ordinality remained relatively constant. Therefore to display results we randomly selected three graphs from each of the scenarios, using the sample function in R (Figure 1). Each of the graphs from the 100 iteractions for each scenario can be found on the author's website. Figure 1: fSVA improves prediction accuracy of simulated datasets. We created simulated datasets (consisting of a database and new samples) using model 2 and tested the prediction accuracy of these using R. For each simulated data set we performed either exact fSVA correction, fast fSVA correction, SVA correction on the database only, or no correction. We performed 100 iterations on each simulation scenario described in table 4. Using the sample function in the R package, we randomly chose three iterations from each of the three scenarios to display here.
We found that fSVA improved the prediction accuracy in all of our simulations ( Figure 1). Interestingly, exact fSVA generally outperformed fast fSVA at all of the correlation levels except the highest correlation levels. However both fSVA methods out-performed our control of performing SVA on the database alone. Additionally any method of batch-correction generally outperformed no batch correction whatsoever.
When the batch and outcome were not correlated, we saw ambiguous performance from using fSVA. This is not unexpected since it has been shown that in scenarios with no confounding between batch and outcome, batch has a minimal effect on prediction accuracy [22]. When databases had extreme confounding between batch and outcome

Results from Microarray Studies
We examined the effect that fsva had on several microarray studies, obtained from the Gene Expression Omnibus (GEO) website [6]. All except three of the studies were preprocessed/standardized as described previously. Three of the studies (GSE2034, GSE2603, GSE2990) were obtained from GEO and fRMA-normalized.
Each of the studies was randomly divided into equally-sized "database" and "new sample" subsets. We SVA-corrected the database subset, and then built a predictive model (PAM) on that corrected data. We then performed fSVA correction on the new samples. After performing fSVA correction, we measured the prediction accuracy of the model built on the database. This process was iterated 100 times for each study to obtain confidence intervals. This method is virtually identical to the simulation described above.
Results from this process can be found below (  Table 2: fSVA improves prediction accuracy in 5 of the 9 studies examined. The remaining 4 studies showed indeterminate results since the 95% confidence intervals overlapped zero. In order to find the prediction accuracy results, each of the studies was randomly divided into "database samples" and "new samples". fSVA-correction was then performed as described above. We then built a predictive model (PAM) on the database and tested the prediction accuracy on the new samples.

Conclusions
Batch effects have been recognized as a crucial hurdle for population genomics experiments [18,22]. They have also been recognized as a critical hurdle in developing genomic signatures for personalized medicine [21]. Here we have introduced the first batch correction method specifically developed for prediction problems. Our approach borrows strength from a training set to infer and remove batch effects in individual clinical samples.
We have demonstrated the power of our approach in both simulated and real gene expression microarray data. However, our approach depends on similarity between the training set and the test samples, both in terms of the genes affected and the estimated coefficients. In small training sets, these assumptions may be violated. Similarly, training sets that show near perfect correlation between batch variables and biological classes represent an extreme case that can not be directly corrected using fSVA. An interesting avenue for future research is the use of publicly available microarray data to build increasingly large training databases for batch removal.
The methods we have developed here are available as part of the sva Bioconductor package [19].