It's not a classical interaction, which is like you say, a cross between two different groups (like treated vs control for cells A and B). But algebraically (A - B) - (C - D) will give you the genes where the differences between A and B are different from the differences between C and D, which it seems is what you are after?

As Gorden says below, doing this test in limma or edgeR is routine, because you have to directly specify the design matrix and contrast. I think you can also do it using DESeq2, by doing something like

design <- model.matrix(~ 0 + gl(4,3))
colnames(design) <- LETTERS[1:4]
design(dds) <- design
dds <- DESeq(dds)
results(dds, contrast = c(1,-1,-1,1))

If I understand your aim correctly, you want something like

pheno.ext = data.frame( AB.vs.CD = as.numeric(Experiment_Groups %in% c(A, B)),
                       AC.vs.BD = as.numeric(Experiment_Groups %in% c(A, C)))

Then run DESeq2 or limma on your data, the pheno.ext data frame as colData (you may need to assign rownames to make DESeq2 happy) and with design ~AB.vs.CD * AC.vs.BD. Then take the results for the coefficient AB.vs.CD:AC.vs.BD. I hope I didn't mess up in this, definitely check that the results make sense (significant genes with a positive coefficient should behave as you expect, higher difference of A vs. B than C vs. D).