Batch effect

Merging two gene expression studies, ComBat

[math]\displaystyle{ \begin{align} Y_{ijg} = \alpha_g + X \beta_g + \gamma_{ig} + \delta_{ig} \epsilon_{ijg} \end{align} }[/math] where X consists of covariates of scientific interests, while [math]\displaystyle{ \gamma_{ig} }[/math] and [math]\displaystyle{ \delta_{ig} }[/math] characterize the additive and multiplicative batch effects of batch i for gene g.

The batch corrected data is [math]\displaystyle{ \begin{align} \frac{Y_{ijg} - \hat{\alpha_g} - X \hat{\beta_g} - \hat{\gamma_{ig}}}{\hat{\delta_{ig}}} + \hat{\alpha_g} + X \hat{\beta_g} \end{align} }[/math]

• Statistics for Genomic Data Science (Coursera) and https://github.com/jtleek/genstats
• Some possible batch variables: operators, runs, machines, library kits, laboratories.
• Alternative empirical Bayes models for adjusting for batch effects in genomic studies Zhang et al. BMC Bioinformatics 2018. The R package is BatchQC from Bioconductor.
• sva::ComBat() function in sva package from Bioconductor.
  • The original paper Johnson 2007 is number 2 of highly cited articles. Figure 1 shows the EB adjustment has the advantage of being robust to outliers in small sample sizes (batch sizes are small).
  • It can remove both known batch effects and other potential latent sources of variation.
  • The tutorial includes information on (1) how to estimate the number of latent sources of variation, (2) how to apply the sva package to estimate latent variables such as batch effects, (3) how to directly remove known batch effects using the ComBat function, (4) how to perform differential expression analysis using surrogate variables either directly or with thelimma package, and (4) how to apply “frozen” sva to improve prediction and clustering.
  • Figure 1 shows 3 heatmaps. Each contains column annotation including Time, Treatment and Batch variables. a) No adjustment, b) standardize each gene within each batch (implemented in dChip software), c) EB batch adjustment. Note that there is no strong evidence of batch effects after adjustment in heat maps (b)–(c).
  • Figure S1 shows the principal component analysis (PCA) before and after batch effect correction for training and validation datasets from another paper
  • Tutorial example to remove the batch effect

BiocManager::install("sva")
library(sva)
library(bladderbatch)
data(bladderdata)
pheno = pData(bladderEset)
edata = exprs(bladderEset)
batch = pheno$batch
table(pheno$cancer)
# Biopsy Cancer Normal 
#      9     40      8 
table(batch)
# batch
#  1  2  3  4  5 
# 11 18  4  5 19 

modcombat = model.matrix(~1, data=pheno)
combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, 
                      prior.plots=FALSE)
# This returns an expression matrix, with the same dimensions 
# as your original dataset (genes x samples).
# mod: Model matrix for outcome of interest and other covariates besides batch
# By default, it performs parametric empirical Bayesian adjustments. 
# If you would like to use nonparametric empirical Bayesian adjustments, 
# use the par.prior=FALSE option (this will take longer). 

combat_edata = ComBat(dat=edata, batch=batch, ref.batch=1)

• ref.batch for reference-based batch adjustment. mean.only option if there is no need to adjust the variancec. Check out paper Alternative empirical Bayes models for adjusting for batch effects in genomic studies Zhang 2018. Figure 4 shows reference-based ComBat can clearly show the pathway activated samples in Batch 1 samples and show the true data pattern in Batch 2 samples from the simulated study (vs the original ComBat approach failed for both cases). In Figure 5 when we cluster genes using K-means, referenced-based Combat can better identify the role of DE or control genes (compared to the original ComBat method).
• Merging two gene-expression studies via cross-platform normalization by Shabalin et al, Bioinformatics 2008. This method (called Cross-Platform Normalization/XPN)was used by Ternès Biometrical Journal 2017.
• Batch effect removal methods for microarray gene expression data integration: a survey by Lazar et al, Bioinformatics 2012. The R package is inSilicoMerging which has been removed from Bioconductor 3.4.
• Question: Combine hgu133a&b and hgu133plus2. Adjusting batch effects in microarray expression data using empirical Bayes methods
• removeBatchEffect() from limma package
• Batch effects and GC content of NGS by Michael Love
• 困扰的batch effect
• Some note by Mikhail Dozmorov

MultiBaC- Multiomic Batch effect Correction

MultiBaC

Combat or limma?

Batch effects : ComBat or removebatcheffects (limma package) ?