My issue is that my samples have been sequenced in 4 batches so I would like to remove batch effects.
Applying ComBat to matrix M (see section 4.7.4 of the manual) indeed seems to work n the MDS plots.
How can I incorporate this change though, to my DGEList object before I proceed to differential methylation analysis is unclear to me.
Can I just apply ComBat to my DGEList counts matrix?
Should I normalize my data any other way apart from what is described in section 4.7.3 of the manual?
ComBat can't be run on a DGEList, but it is not necessary to do so.
The best method to correct for batch effects is to include batch as a factor in the edgeR linear model. Then edgeR will correct for any batch effects as part of the differential methylation analysis.
In my opinion, it is always better to adjust for batch effects as part of the linear model (see Section 3.4.3 of the edgeR User's Guide) rather than as a separate step. Running batch correction as a separate step tends to change the data in risky ways and can very often result in spurious significant results.