svaseq: removing batch effects and other unwanted noise from sequencing data

1Batch effects and other technological artifacts introduce spurious correlation, create bias and add variability to the results of genomic experiments. The basic problem is that batch effects introduce a new source of signal into the data that can be confused with the signal an analyst is looking for. This signal is consistent across transcripts, exons or genes and so may lead to gross errors in the calculation of statistical significance, estimates of effect sizes or other statistical measures. These types of noise also prevent analysts from appropriately modeling biological variation and group-specific changes in gene expression. Unfortunately we rarely know all of the potential artifacts in most high-throughput experiments. In some cases, it is possible to rely on the date the samples were processed as a surrogate for unmeasured artifacts and correct for them to get statistically accurate results. However, each new technology may suffer from different artifacts and it may take time for the community to discover which variables must be measured and included in an analysis.

2In 2007 we introduced surrogate variable analysis (sva) as a conceptual approach to statistical modeling of genomic data when artifacts are unknown or unmeasured (Figure 1) and subsequently improved the estimation algorithm (Figure 2). We proposed modeling the data as a combination of known variables of interest, known adjustment variables and unknown and unmeasured artifacts. A simple version of this model might relate gene expression for gene i on sample j (gij) to the phenotype for that sample (yj), the known batch variable for that sample (aj) and an unknown artifact on that sample (uj):

3If only a single gene is measured, it is difficult to estimate the unknown artifact (uj) from the data directly, since all the coefficients (b, c, d) are unknown. But we noticed that if many genes are measured, it is possible that for some genes the coefficients for bi1 and ci may be equal to zero. For these specific genes the model reduces to: Our next insight was that even though di and uj are unknown, you do not need to know either of them exactly to get correct statistical inference for the parameters for the phenotype variable bi1, you just have to know their linear combination di × uj. We showed that if you collect the data for all genes where there is no effect of phenotype or known batch (bi1 = 0 and ci = 0) and subtracted the mean of each gene to remove the baseline effect, the matrix form of the model is: where mu is the number of genes where b1i = 0 and ci = 0 and n is the sample size. We then pointed out that you could apply a matrix decomposition like the singular value decomposition or principal components analysis to this subset of the data to estimate surrogates for the unknown artifacts u. Later we also showed that if the number of genes that are not affected by phenotype or known batch (b1i = 0 and c = 0) but are affected by unknown artifacts (di ? 0) is large enough, you can obtain consistent estimates of a linear transformation of ui. You can then include the estimate of batch effects in downstream models to remove artifacts, dependence between genes and improve statistical inference. These insights inspired the general form of the sva approach shown in Figure 1. The difference between the two-step and iteratively re-weighted sva algorithms is the approach to estimating the genes affected by unknown artifacts (Step 1).

4In this paper I discuss two new extensions of the sva approach for sequencing data. The first extension deals with Step 1 of the sva approach (Figure 2). The idea is to use external information to estimate the genes that are likely to be affected by unknown artifacts. The original sva algorithm attempted to identify genes affected by artifacts directly from the data itself, but sometimes there is external information about which genes are unlikely to be differential expressed. This external information could be control probes or estimates of batch-related probes from previous studies. Supervised sva uses this information in Step 1 of the sva approach, reducing computational time and reliance on estimation procedures for identifying the right genes to use for artifact estimation. The second idea, svaseq, is based on performing an appropriate transformation of the count or Fragments Per Kilobase Of Exon Per Million Fragments Mapped (FPKM) data during Steps 1 and 2 of the sva approach. Here I focus on the moderated log transform, which has been widely adopted both for the analysis of sequence count data and FPKM estimates. I then perform a thorough and reproducible comparison between the standard methods for removing batch effects from sequencing data. I demonstrate that svaseq and supervised sva perform comparably to existing approaches for removing batch effects from sequencing data.

5The general form of the simple model in Equation (1) is or in the matrix form is where the function f has been applied component wise to each element of G and there may be multiple phenotypes, artifacts or unknown artifacts and I have dropped the explicit conditioning for ease of notation. The general sva algorithm then proceeds in the following four concrete steps.

6In Step 5, the covariates may be included in a standard linear regression modeling analysis on an appropriately transformed scale or the covariates can be directly used in software that models counts with generalized linear models (GLMs) including edgeR and DESeq. They can be directly included in these models as they are estimated on the same scale as standard link functions for GLMs.

7The general form of the sva algorithm relies on the idea that there is a subset of genes, probes or transcripts that are affected by unknown batch effects or other artifacts but are not affected by the biological relationship of interest. This is the main way that a sva approach is distinguished from a standard principal component regression. Standard principal component-based regression methods, such as EIGENSTRAT, are sufficient when the number of probes, genes or features expected to show a signal is small. Then the principal components will be consistent estimates of a linear transformation of the artifacts and not the phenotype/outcome of interest.A number of extensions and variations on the sva algorithm have been introduced. For Step 1 in the sva algorithm, identifying probes only associated with unmeasured artifacts, it has been proposed to use control probes. For Step 2 of the sva algorithm, estimating latent factors only associated with unmeasured artifacts, it has been proposed to use factor analysis, independent components analysis and principal components analysis. Another extension of the surrogate variable approach in Step 2 has been to model known sources of technical or biological covariation between the measurements for probes, for example in eQTL studies.

8Supervised sva (ssva) sets λi = 1 for all negative controls and λi = 0 for all other genes in Step 2 of the sva algorithm. The assumption is that control probes will capture all of the variation due to unknown artifacts and none of the variation due to the phenotype. Control probes may miss biological artifacts. For example, we showed that trans-eQTL that are associated with multiple gene expression levels may act like an artifact when measuring the association between gene expression and phenotype. These artifacts may be missed by the ssva approach. However ssva is particularly useful for unfortunate experimental designs where the phenotype variable and unknown artifacts are highly correlated, making empirical estimates unstable.

9The second extension involves the choice of function f in (1). In our original work, we used the identity function for data measured on an approximately symmetric and continuous scale. For sequencing data, which are often represented as counts, a more suitable model may involve the use of a moderated log function. For example in Step 1 of the algorithm we may first transform the gene expression measurements by applying the function log(gij + c) for a small positive constant. In the analyses that follow I will set c = 1. The choice of the moderating constant is an important one and is beyond the scope of this manuscript. Intuitively a choice of c = 0 corresponds to no moderation and as c increases you decrease the variation in the data. After performing Steps 1–5 of the sva estimation algorithm, the estimated covariates are included in downstream models as adjustment variables. For the analyses that follow, I will use the limma package with the voom method for differential expression analysis. The voom method is an approach for estimating the mean–variance relationship when performing differential expression analysis on sequencing experiments.

10Supervised svaseq proceeds by applying the transformation log(gij + c) to the gene expression count data in Step 1 and setting λi = 1 for all negative controls and λi = 0 for all other genes in Step 2 of the sva algorithm.

11I use data from Zebrafish sensory neurons with three control samples and three gallein treated samples as the comparison groups. These data are available as part of the zebrafishRNASeq Bioconductor package. I loaded the data and filtered as described in the removing unwanted variation in sequencing data (RUVSeq) package. Then I estimated batch effects using supervised and unsupervised sva for sequencing, principal components analysis, RUV with control probes, RUV with empirical controls and residual RUV. I compared the model estimates and I compared differential expression analysis results when each of the different batch effect estimates was included in the model in place of the study variable.

12ReCount is a database of pre-processed RNA-sequencing data, processed to be comparable across samples. In this analysis, I downloaded pre-counted RNA-sequencing datasets measuring gene expression in two separate Hapmap populations. For my analysis, I downloaded the count data from ReCount and downloaded the pedigree information from the Hapmap website. I then performed differential expression analyses looking for differences in expression between males and females and estimated unknown latent structure. I calculated estimates of batch effects using unsupervised sva for sequencing, principal components analysis, RUV with empirical control probes, and RUV on residuals. I compared the estimates to the variable indicating whether the data came from the Pickrell or Montgomery study. I compared two scenarios, one where the sex and study variables were balanced and one where they were imbalanced. I also compared the differential expression analysis results when each of the different batch effect estimates was included in the model in place of the study variable.

13I downloaded the processed GEUVADIS Ballgown object from:https://github.com/alyssafrazee/ballgown_codeI then subset the data to only the non-duplicated samples and performed a differential expression analysis comparing populations. I calculated estimates of batch effects using unsupervised sva for sequencing, principal components analysis, RUV with empirical control probes and RUV on residuals. I compared the estimates to the variable indicating which lab the sequencing was performed in. I also compared the differential expression analysis results when each of the different batch effect estimates was included in the model in place of the laboratory variable.

14I simulated data from a negative binomial model for count based RNA-sequencing data (Figure 3). For complete details see the simulated data R markdown document and accompanying HTML file.I estimated the model parameters from the Zebrafish data described above. I simulated two scenarios, one where the group and batch variable were not correlated and one where they were correlated. Here we consider both batch effects that are correlated with the outcome and batch effects that are orthogonal. This is a critical distinction as unsupervised methods that estimate batch effects directly from the data will often perform worse when batch and outcome are correlated unless this relationship is explicitly modeled. Data were simulated with the Polyester R package. Then I estimated batch effects using supervised and unsupervised sva for sequencing, principal components analysis, RUV with control probes, RUV with empirical controls and residual RUV. I compared the model estimates to the true simulated batch variable and I also compared differential expression analysis results when each of the different batch effect estimates was included in the model in place of the study variable.

15ssva and svaseq are currently implemented in the sva Bioconductor package version 3.11.2 or greater (http://www.bioconductor.org/packages/devel/bioc/html/sva.html). All data and code used to perform this analysis are available as R markdown files available from: https://github.com/jtleek/svaseq. You can view the individual analyses as webpages at: Zebrafish analysis: http://jtleek.com/svaseq/zebrafish.htmlReCount analysis: http://jtleek.com/svaseq/recount.htmlGEUVADIS analysis: http://jtleek.com/svaseq/geuvadis.htmlSimulated data analysis: http://jtleek.com/svaseq/simulateData.html

16Zebrafish analysis: http://jtleek.com/svaseq/zebrafish.html

17ReCount analysis: http://jtleek.com/svaseq/recount.html

18GEUVADIS analysis: http://jtleek.com/svaseq/geuvadis.html

19Simulated data analysis: http://jtleek.com/svaseq/simulateData.html

20I estimated simulation parameters from the Zebrafish data as described in the methods. I then performed several checks to confirm that (i) the data generated by the simulated model recapitulated the qualitative behavior of the data used to estimate the model parameters (Figure 4), (ii) that data generated without signal did not show statistically significant results, (iii) data could be simulated with differential expression signal and (iv) that data with batch effects displayed the expected conservative bias of P-values (See supplementary analysis files (http://jtleek.com/svaseq/zebrafish.html).

21In the Zebrafish study there was a single simulated batch effect. I simulated two scenarios, in the first scenario the batch effect and the group effect had low correlation. As expected, all methods that aim to estimate batch effects while taking into account multiple sources of signal (svaseq and RUV methods) produce estimates that are highly correlated with the simulated batch effect. The estimate of batch based on principal components is biased, because the principal component is estimating a linear combination of the group and batch variable (Figure 5). In this scenario all the P-value distributions are approximately correct, with the exception of the analysis using principal component based estimates of batch. This is because the principal components do not estimate accurate versions of the batch effect and bias the statistical significance calculation (see: http://jtleek.com/svaseq/simulateData.html for plots).

22I next fit models relating the gene expression counts to the simulated group phenotype (p) and batch effect estimates using the following model: I accounted for the potential relationship between mean and variance using the voom method. I then estimated how concordant the differential expression results were with the results we obtained when we fit model (6) using the true simulated batch variable using concordance at the top plots (CAT plots, Figure 6). For each ranking, these plots show the fraction of results that are the same between the analysis using the true batch variable and the batch variable estimated with different methods. Supervised and unsupervised sva for sequencing, RUV with control probes, RUV with empirical controls and residual RUV all outperform not adjusting for batch effects while principal components analysis performs worse than no adjustment. The reason is that the principal component estimate is correlated with group and absorbs some of the signal due to that variable.

23We performed an identical analysis where the batch was now correlated with the group variable. Qualitatively similar results hold in this second simulated scenario with one exception. The empirical RUV methods attempt to define control probes by identifying genes that do not show differential expression with respect to batch. But when batch and group are correlated, this may also through away genes that show signal with respect to the group variable. Similar the residual RUV approach estimates the batch variable after taking the residuals from the model fit of the counts on the group variable. However, when batch and group are correlated, this again removes batch signal and leads to slightly lower performance of the RUV approaches. Unsupervised svaseq does not use the control probes but avoids some of these difficulties by iteratively identifying probes associated with group but not associated with batch (Figure 7). The P-value histograms here show a strong difference between supervised and unsupervised approaches. The unsupervised approaches attempt to estimate the artifacts, but they are correlated with the group. Since the estimates are off, the statistical significance calculations are not correct (see: http://jtleek.com/svaseq/simulateData.html for plots). But the supervised methods correct the statistical significance calculations accurately.

24Next I performed an analysis on the Zebrafish data as described in the methods section. Here, the batch variable is not known, but we do have negative control probes which can be used to estimate the batch effects. When comparing the batch estimates, I noted that the supervised sva estimates and the RUV control probes estimates were perfectly correlated (R2 = 1) and that they produced identical differential expression results (Figure 8a). The unsupervised sva and principal components approaches are most similar to the supervised estimates from SVA or RUV for the Zebrafish data.

25For the ReCount data I generated an artificial batch effect by combining the data from two different studies of gene expression in two different populations. I used sex as the outcome variable in the analysis and then estimated batch effects using the same set of proposed approaches. I next fit models relating the gene expression counts (g) to sex variable (p, phenotype variable representing sex) and batch effect estimates (, representing the study the samples came from) using the following model: I accounted for the potential relationship between mean and variance using the voom method. I then estimated how concordant the differential expression results were with the results we obtained when we fit model 6 using the true simulated batch variable using concordance at the top plots (Figure 8b).

26In the original data, the batch effect and the group variable are nearly perfectly orthogonal. In this situation the empirical and residual RUV approaches produce estimates of the batch variable more highly correlated than the unsupervised svaseq approach (see http://jtleek.com/svaseq/recount.html) and produce correspondingly more similar differential expression results to using the true study variable as an adjustment in the differential expression analysis (Figure 9a). However, I next re-sampled the data to mimic a scenario where the group and batch variable showed modest correlation (r2 = 0.33). In this scenario the unsupervised sva and principal components analysis approaches outperform the empirical control RUV approach. The residual RUV approach performs worse than no adjustment for study, because signal due to the batch variable was removed when the residuals from the model relating sex to phenotype was calculated (Figure 9b).

27The Ballgown R package https://github.com/alyssafrazee/ballgown can be used to analyze abundance data from assembled transcriptome data from Cufflinks. I loaded data from the GEUVADIS project that we recently processed using Cufflinks and Ballgown. I selected only the non-duplicated samples and performed a differential expression analysis comparing different populations. I then compared the estimated batch effects using the various approaches to the known lab where the samples were processed, one of the variables that showed the highest association with assembled transcript levels.I assessed concordance between the batch effect estimates and the lab variable by fitting the model: and then performed an ANOVA to compare the model including lab to the null model of no association with lab. The unsupervised sva and principal components estimates showed significantly higher F-statistics for concordance (482 and 456, respectively) compared to the RUV approach (106 and 109 for RUV residual and empirical, respectively). I next fit models relating the gene expression counts (g) to the population phenotype (p, phenotype variable representing population) and batch effect estimates (, representing the study the samples came from) using the following model: I compared the results to the differential expression model where I included the known lab variable as an adjustment in place of. The svaseq and principal components adjusted analyses showed greater concordance with the lab adjusted analysis, as expected since the batch estimates were more highly correlated with this known variable (Figure 10).

28Here I have described the general sva framework and I have introduced two extensions of the sva approach. The first takes advantage of known control probes to simplify the sva algorithm and the second addresses the distribution of count and FPKM data typically observed in sequencing experiments. The question of whether to use FPKM or count based approaches for the analysis of RNA-sequencing data is beyond the scope of this paper. However, I have demonstrated in this paper that regardless of the choice for measurement summary, svaseq can be applied to remove batch effects.I have shown that sva-based approaches perform comparably to other batch effect estimation procedures for sequencing when the group and unknown batch variables are uncorrelated and outperform other approaches when the batch and group variable are correlated. These extensions are currently available from the devel branch of the sva software http://bioconductor.org/packages/devel/bioc/html/sva.html and all analyses are fully reproducible and available as R markdown documents from https://github.com/jtleek/svaseq.