Dear all,

I analyzed a batch of samples using edgeR in order to isolate differentially expressed genes from a bacteria strain. Two experimental conditions were varying (growing phase and temperature), and I had at least 3 replicates for every set of conditions. We are mostly interested in temperature effect but I also extracted DE genes between different growing phases (with blocking of temperature condition).

Everything went fine until I received a second batch of data with additional temperature values: one control condition (temp1, same as in batch01) for consistency check, and a new temperature value. QC look fine (reads quality, alignment, etc...), and a quick look with a genome browser showed no sign of excessive noise (see picture for example). However, samples from second batch were sequenced much deeper (bottom, ~25 million reads) than ones from first batch (top, ~7 million reads).

In order to check for consistency between batches, I generated an MDS plot from normalized cpm counts (see figure).

From this picture, we can see that samples for control condition ('Temp1', common between 2 batches) do not cluster together between batches (small VS large points) for 'G2'. I suspect a pronounced batch effect.

If you confirm that this looks like a batch-effect problem, what would be the best strategy for correcting it (if possible) before EdgeR analysis ?

Is it a problem to have only a single temperature in common between batches ?

I read about limma's "removeBatchEffect()" function but I would like confirmation that I can use it in this context. Actually, I suspect that it would forbid any use of my edgeR pipeline because of the transformation from counts to continuous data.

I also found references to Voom batch effect modeling, would it be more relevant than any EdgeR approach in this context ?

Thank you very much for your help.

R.

This does look like a batch effect. Fortunately, it seems you have groups that are present in both batches (i.e., Temp1 G1 and G2). This means that the experimental design is not completely confounded. You can include a batch factor in your design matrix, and edgeR will automatically consider the batch effect when performing comparisons. This is possible as the batch effect is effectively estimated from the inter-batch differences within the Temp1 G1/G2 groups. There are a number of examples of blocking factors in the edgeR user's guide, see Section 3.4 for details.

Whatever you do, do not try to remove the batch effects from the raw data prior to running edgeR. Functions like removeBatchEffect are intended for removing the batch effect prior to using model-free methods (e.g., clustering, dimensionality reduction). It is much better to model the batch effect explicitly by including a blocking factor in the design matrix, as this ensures that the uncertainty of the batch effect estimate is properly considered during hypothesis testing.