Hi,

I have been using Combat to correct the batch effect in 450k data (~10 batches). Recently, I read an old reply from Dr. Peter Langfelder where he mentioned that  "ComBat should NOT be used before running association testing (lmFit); association testing should be run with batch as a covariate on original data."

I have read the paper from Combat authors where they mentioned that Combat performs better than SVD when sample size is small and comparably similar for large sample size.

Now, my question is since my sample size is large and I am using limma for calculating differential methylation, should I be adjusting batch directly in limma as a covariate (or use removeBatchEffect(), not sure about this, as this will again be like removing batch separately) 

batch <- pheno$batch
BSC <- pheno$BSC_batch
group <- factor(targets$status,levels=c("Control","Case"))
design <- model.matrix(~targets$Age+batch+BSC+group)
fit <- lmFit(Mval, design)

or should I continue with Combat for removing batch and limma for adjusting other biological covariates?

MvalC <- ComBat(Mval, batch=batch, mod=NULL, par.prior = TRUE, prior.plots = FALSE)
modcombat <- model.matrix(~target$Age, data=pheno)
MvalC1 <- ComBat(MvalC, batch=BSC, mod=NULL, par.prior = TRUE, prior.plots = FALSE)
design <- model.matrix(~targets$Age+group)
fit <- lmFit(MvalC1, design)

(I corrected for Chip batch and BSC batch separately, as they were confounded and limma was showing error.)

My recommendation is to try using ComBat and then also including a batch covariate in your model. This accounts for any residual df problems (mentioned above) and allows you to use the extra benefits of ComBat that are not available by just using a blocking factor (i.e. adjusting for differences in batch variance). 

To elaborate further: if there are no variance differences between the batches, adding a blocking factor is just fine--there will be no difference in using ComBat or a blocking factor. However, this is not common--in most cases the batch effect impacts both the means and variances across batches. The greater the differences in variance between batches, the greater benefit you will see from using ComBat+Limma(plus batch covariate) versus Limma (plus batch) alone.  

My colleagues and I spent considerable time looking into problems with attempting to remove the batch effect, which resulted in the article:
http://biostatistics.oxfordjournals.org/content/17/1/29

I do not disagree with the other answers from Evan and Aaron. I will explain using somewhat different words why you should avoid your second alternative.

ComBat, removeBatchEffect and similar methods can not truly remove the batch effect, nor can they create a batch effect free data set. For that to happen, the true batch effect must be known, which rarely is the case. Instead, as in your example, it is estimated from the data and it is this estimated batch effect that is removed from the data. Thus the difference between the true batch effect and the estimate, i.e. the estimation error, will be present in your adjusted data. This estimation error is in itself a batch effect. So this procedure swaps one batch effect with another, one problem with another quite similar problem. It can be hard to tell to which degree the second batch effect is detrimental, but be aware that it will be there even if you start out with no true batch effect.

So your second analysis alternative is problematic because when you use ComBat you ensure your data has a batch effect (the estimation error). When this batch effect is ignored in your subsequent limma analysis, unreliable p-values may be the result (depending on other circumstances like experiment design)

Your first alternative is more a safe textbook-like way of doing it.