Batch effect

Merging two gene expression studies, ComBat

• Statistics for Genomic Data Science (Coursera) and https://github.com/jtleek/genstats
• Some possible batch variables: operators, runs, machines, library kits, laboratories.
• sva::ComBat() function in sva package from Bioconductor. [math]\displaystyle{ \begin{align} Y_{ijg} = \alpha_g + X \beta_g + \gamma_{ig} + \delta_{ig} \epsilon_{ijg} \end{align} }[/math] where [math]\displaystyle{ X=X_{ij} }[/math] consists of covariates (eg biological) of scientific interests (e.g. Pathway activation levels in Zhang's 2018 simulation example), 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 error terms, [math]\displaystyle{ \epsilon_{ijg} }[/math], are assumed to follow a normal distribution with expected value of zero and variance [math]\displaystyle{ \sigma^2_𝑔 }[/math]. 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]
• Alternative empirical Bayes models for adjusting for batch effects in genomic studies Zhang et al. BMC Bioinformatics 2018. The R package is sva and BatchQC from Bioconductor.
  • Reference batch adjustment: [math]\displaystyle{ \begin{align} Y_{ijg} = \alpha_{rg} + X \beta_{rg} + \gamma_{rig} + \delta_{rig} \epsilon_{ijg} \end{align} }[/math] where [math]\displaystyle{ \alpha_{rg} }[/math] is the average gene expression in the chosen reference batch (r). Furthermore, [math]\displaystyle{ \gamma_{rig} }[/math] and [math]\displaystyle{ \delta_{rig} }[/math] represent the additive and multiplicative batch differences between the reference batch and batch i for gene g. The error terms, [math]\displaystyle{ \epsilon_{ijg} }[/math], are assumed to follow a normal distribution with expected value of zero and a reference batch variance [math]\displaystyle{ \sigma^2_{𝑟𝑔} }[/math].
  • Mean-only adjustment for batch effects: [math]\displaystyle{ \begin{align} Y_{ijg} = \alpha_{g} + X \beta_{g} + \gamma_{ig} + \epsilon_{ijg} \end{align} }[/math]
• svg vignette 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). In addition to the github reposition for the simulation R code, BatchQC::rnaseq_sim() can also do that.
  • Paper source code, Simulations on reference batch & modsel_sensitivity_uponly.R.
• 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
• Figure S1 shows the principal component analysis (PCA) before and after batch effect correction for training and validation datasets from another paper
• Batch effects and GC content of NGS by Michael Love
• 困扰的batch effect
• Some note by Mikhail Dozmorov

ComBat-Seq

• sva package vignette

  # Example 1
  count_matrix <- matrix(rnbinom(400, size=10, prob=0.1), nrow=50, ncol=8)
  batch <- c(rep(1, 4), rep(2, 4))
  adjusted <- ComBat_seq(count_matrix, batch=batch, group=NULL)
  
  # Example 2 -  one biological variable
  group <- rep(c(0,1), 4)
  adjusted_counts <- ComBat_seq(count_matrix, batch=batch, group=group)
  
  # Example 3 -  multiple biological variables
  cov1 <- rep(c(0,1), 4)
  cov2 <- c(0,0,1,1,0,0,1,1)
  covar_mat <- cbind(cov1, cov2)
  adjusted_counts <- ComBat_seq(count_matrix, batch=batch, group=NULL, covar_mod=covar_mat)
  

• ? ComBat_seq It currently (sva version 3.44.0 as of 2022/6) does not have the ref.batch option? See this post Can you add the parameter of ref.batch in ComBat-Seq? #51. One alternative is to normalize the count data of the two batches separately first (for example, using DEseq2) and correct the batch effect then using combat.
• ComBat-seq: batch effect adjustment for RNA-seq count data 2020
• Introduction to Bioconductor SVA and ComBat-Seq in R
• Using ComBat-seq on transcript counts.
  • You need to use the RSEM expected counts (values in the expected_counts column are still not integers).
  • There is no need to round them to exact integers. You absolutely cannot use TPM or FPKM.
  • An example of RSEM output here.
• How to use combat in order to remove batch effects? DESeq() + vst() + limma::removeBatchEffect()
• 数据分析：RNA-seq数据的批次校正方法 It includes complete R code for running ComBat_seq() including ComBat_seq(), limma+removeBatchEffect, PCA, VOOM+SNM, DESeq2.
• Evaluation of critical data processing steps for reliable prediction of gene co-expression from large collections of RNA-seq data Vandenbon 2022

DESeq2

• Batch effects vs biological variables
• Batch-correction on tximport data from salmon. Include the batch effect in the design of DESeq2. No need to use an external tool like ComBat_seq(). Search the keyword batch in the Vignette.
• Batch correction for RNA-seq didn't work with ComBat-seq tool

limma::removeBatchEffect()

• removeBatchEffect() from limma package
• Batch effects : ComBat or removebatcheffects (limma package) ?
• Tutorial:Basic normalization, batch correction and visualization of RNA-seq data

MultiBaC- Multiomic Batch effect Correction

MultiBaC