I am using the R (using EdgeR) for the RNA Seq analysis, I had few batch effect samples like Control vs treatment. Could anyone tell me the best way to remove the batch effects.

I have looked into the cluster dendrogram and MDS post (using EdgeR), the Control and Treatments are not grouping together, and their is differences in the different batch within Control and Treatment.

I tried this syntax to remove the batch effects for the rest of the analysis to get the DE Genelist. I feel it didn't able to remove the batch effects from the rest of the analysis.

# define the experiment
Batch <- factor(c(1,2,3,1,2,3))
Treatment <- factor(c(0,0,0,1,1,1))
design_noI <- model.matrix(~ 0 + Batch + Treatment)

# define our DGEList
San <- DGEList(counts=x,genes=gene_info)

A <- aveLogCPM(San)
y2 <- San[A>1,]
y2 <- calcNormFactors(y2)
logCPM <- cpm(y2, log=TRUE, prior.count=5)
logCPMR <- removeBatchEffect(y2, Clutch = Clutch)

I can see the difference in MDS plot on removing the batch effects but I couldn't able to carry the same data flow into my rest of analysis.

# filter the lowly expresed genes and normalize the data
keep <- rowSums(cpm(San/logCPMR)>1) >= 3
keep2 <- removeBatchEffect(keep, Batch)
y_group <- San/logCPMR[keep2, , keep2.lib.size=FALSE]
y_group <- calcNormFactors(San/logCPMR)

I have mention the San/logCPMR, as San take as data without remove the batch effects, whereas I think logCPMR as batch effected removed data frame.

Any suggestion on this more useful to me.

Please let me know, any other information needed to deal with this batch effects.

  • $\begingroup$ I am not sure if I understand exactly what you're asking (and what you are doing in your code). You seem to deviate from edgeR manual for paired or batch effect analysis. Your first 3 lines of code seem okay to me, but then you should fit your data to the model, including your design with batch effect. See 3.4 of edgeR manual. The removeBatchEffect results are only meant for clustering or visualization not the statistical analysis. $\endgroup$
    – benn
    Commented Jul 25, 2018 at 20:48
  • $\begingroup$ Thanks b.nota, My point is ...when their is batch effects, shown in MDS plot or Cluster Dendrogram, re-visualise it by removing batch effects in MDS plot, But could we take the same batch effect removed data for further data analysis to get the DE genelist. As if we don't remove the batch effects from the data, that effects on the final DE genelist, and I get the less number of DE genelist. $\endgroup$
    – San
    Commented Jul 25, 2018 at 20:56
  • $\begingroup$ No you cannot use removeBatchEffect results for DE analysis, follow the guidelines in manual and see answer of Devon. $\endgroup$
    – benn
    Commented Jul 26, 2018 at 10:38

1 Answer 1


The help page for removeBatchEffect() explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the combat() function from the sva package.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.