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.

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
