I got totally confused when looking at edgeR manual somehow:
The are different ways of estimating the dispersions (GLM and a classic approach). Also, one might estimate either CommonDisp or TrendedDisp.
I get totally different results when applying GLM and classic and also when excluding TrendedDisp.
My experiment is quite simple: I have 2 treated and 2 untreated samples and want to find differential expressed genes. 
I also found that I have a batch effect, so that I got batch-values from RUVseq. What would be an applicable method (GLM or classic) to use if I want to intergrate batch effect into the design matrix. Also, when should I use estimateGLMTrendedDisp over estimateGLMCommonDisp?

GLM:

y <- estimateGLMCommonDisp(y, design)

y <- estimateGLMTrendedDisp(y, design)

y <- estimateGLMTagwiseDisp(y, design)
then apply glmFit and glmLRT

classic:

y <- estimateCommonDisp(y)

y <- estimateTagwiseDisp(y)
then apply exact test

EDIT:
also I get different BCV plots when I run either y <- estimateDisp(y, design) or those three commands together

y <- estimateGLMCommonDisp(y, design)

y <- estimateGLMTrendedDisp(y, design)

y <- estimateGLMTagwiseDisp(y, design)
​

​what is quite weird because according to the edger manual they should give the same results.

It's hard to know what's going on without having some diagnostics, but here's some general points. Firstly, you haven't performed trended dispersion estimation in your code with classic edgeR. This will obviously yield different results to the GLM method where you get trended estimates, because the mean-dispersion trend provides a more accurate model of the empirical mean-variance relationship. Secondly, the GLM functions are implemented differently from estimateDisp, so the dispersions that you get out of them will not be exactly equal, even at the best of times. Finally, it's not surprising that exactTest and glmLRT give different results, as they operate off different statistical frameworks. This shouldn't matter for strongly DE genes, but I suspect that it will affect the genes that are borderline significant.

More educated advice would require examination of the BCV plots - if these are substantially different between methods, then the p-values will inevitably be different. The various dispersion estimation functions should (and do) give similar results for well-behaved data. Large differences should only manifest if there's something weird going on. From the information you've supplied, this is not out of the question for your data set, especially if there's a big batch effect in there somewhere.

As for the batch blocking factor, you'll need to add it as an additive term to your design matrix. You've only got 2 residual d.f. in your design, so you can only afford a single blocking term. You can probably add it in by doing something like this:

design <- model.matrix(~batch + grouping)

... and then make sure you drop the last coefficient to test for DE between groups. Obviously, this only works if your batch effect doesn't confound your groupings of interest. If both samples of one group are in one batch, and both samples of the other group are in another batch, then you won't be able to separate your batch effect from genuine DE.

Edit: Actually, Gordon reminded me that estimateTagwiseDisp fits a trend by default, so my first point shouldn't be a problem. The fitting algorithm is different from the GLM methods, though - the former uses a moving average compared to a bin spline for the latter. This will still cause some differences, though I wouldn't have expected anything too dramatic.