Question: Hi there,
after having read the vignette, I am really keen to understand in depth how the Limma functions 'removebatcheffect' actually works.
Looked at the function:

# removeBatchEffect.R

# A refinement would be to empirical Bayes shrink
# the batch effects before subtracting them.

removeBatchEffect <- function(x,batch=NULL,batch2=NULL,covariates=NULL,design=matrix(1,ncol(x),1),...)
# Remove batch effects from matrix of expression data
# Gordon Smyth and Carolyn de Graaf
# Created 1 Aug 2008. Last revised 1 June 2014.
{
if(is.null(batch) && is.null(batch2) && is.null(covariates)) return(as.matrix(x))
if(!is.null(batch)) {
batch <- as.factor(batch)
contrasts(batch) <- contr.sum(levels(batch))
batch <- model.matrix(~batch)[,-1,drop=FALSE]
}
if(!is.null(batch2)) {
batch2 <- as.factor(batch2)
contrasts(batch2) <- contr.sum(levels(batch2))
batch2 <- model.matrix(~batch2)[,-1,drop=FALSE]
}
if(!is.null(covariates)) covariates <- as.matrix(covariates)
X.batch <- cbind(batch,batch2,covariates)
fit <- lmFit(x,cbind(design,X.batch),...)
beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE]
beta[is.na(beta)] <- 0
as.matrix(x) - beta %*% t(X.batch)
}

My understanding is that when you have just **one** **known** source of bias (say, different days in which an experiment was performed), covariates is = null. Thus, the function will be the following:

batch <- as.factor(batch)
contrasts(batch) <- contr.sum(levels(batch))
batch <- model.matrix(~batch)[,-1,drop=FALSE]

once a design matrix has been created..

fit <- lmFit(x,cbind(design,X.batch),...)
beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE]
beta[is.na(beta)] <- 0
as.matrix(x) - beta %*% t(X.batch)

so that's the hard bit for a lay person...
so once you have your 'dds' object (I am using DESeq2), design formula (es. ~condition), Limma's function adds the batch object just created and then it fits a linear model for each gene in the dataset.

1. Is that correct?
2. can anyone explain what is the beta object in this case? I guess the following strings will add the beta variable to the equation?

Sorry if my questions may be very naive for you but I am really trying to understand every single steps and in the meantime trying to study a bit of stats.

Answer: Obviously *limma* should not be used to model the raw or normalised counts from DESeq2 due to the fact that these counts will follow a negative binomial distribution, whereas *limma*'s `removeBatchEffect()` function [and *limma*, generally speaking] will expect that your data is normalised and follows a normal distribution (to be specific: that the residuals after fitting the model follow a normal distribution).

In a practical setting, the following would occur:

1. via whichever means, the user identifies a batch effect in their
data, `batch`
2. user includes `batch` in the design formula when using DESeq2, for
example, `~ condition + batch`
3. when deriving test statistics for the main condition of interest,
`condition`, DESeq2 will, internally, treat `batch` as a covariate
and 'adjust' the test statistics for `condition`, including the
p-values, obviously. We could just as easily derive test statistics for batch, and these would conversely be adjusted for condition. This is how covariate-adjustments work in regression models.
4. user then wishes to use this data for downstream applications; so,
transforms the normalised counts via `rlog(..., blind = FALSE)` or
`vst(..., blind = FALSE)`

-----------
----

Now, we have performed our differential expression analysis using DESeq2 and adjusted for the effect of batch ***in the test statistics output only*** - no modification to our expression data has been made for `batch`, or, indeed, anything else in our design formula. However, for our downstream programs / methods, like clustering, PCA, and other cool stuff, we want to completely eliminate the batch effect directly from our expression data; so, we use:

limma::removeBatchEffect(assay(vst), batch = colData(dds)[,batch'])

Limma, via the code that you have posted above, will estimate the effect of batch for each gene by fitting a model of form `~ batch` [for each gene], and then use these estimates ('beta coefficients') to directly modify the expression levels in your `assay(vst)` object. So, the '*beta*' to which you inquire is a measure of the effect of batch for each gene. In a typical regression scenario, we would derive odds ratios from these beta coefficients - the odds ratio is the exponent of the beta coefficient.

Trust that this helps.

Kevin