Looking at this, the groups structure does seem to be wrong for what you want;
groups = list(NDE = c(1,1,1,1,1,1), DE = c(1,1,2,2,1,1)))
indicates that the two models under consideration are one in which all samples behave identically, and one in which the third and fourth samples behave differently to the first, second, fifth and sixth.
You're also missing a libsizes declaration (use the getLibsizes function, if you don't have an alternative way of filling this slot).
However, you want to find differences between the timepoints, with no differences between individuals. There are a couple of ways to do this in baySeq.
The simplest way is to ignore any patient specific effects and treat each timepoint as a separate replicate group. (This is the correct approach for destructive sampling of samples at different timepoints, which is often the case). Then your code would look like:
cD <- new("countData", data = otus[,2:25], replicates = rep(1:4, each = 6),
groups = list(NDE = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
DE = c(1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1)))
and then proceed as usual with a negative binomial distribution (the default). Alternatively, you could use the 'allModels' function and consensus priors approach to check diverse patterns behaviour over timepoints; see section 3.4 of the preprint http://biorxiv.org/content/early/2014/11/28/011890.
The second way is to try and incorporate patient specific effects into a multinomial-dirichlet distribution. Here the code would look more like what you already have, though the groups slot changes:
cD <- new("countData", data = cdat, replicates = c(1,1,1,1,1,1), groups = list(NDE = c(1,1,1,1,1,1))
libsizes(cD) <- getLibsizes(cD)
densityFunction(cD) <- mdDensity # use multinomial-dirichlet distribution
cD <- getPriors(cD, cl = cl)
cD <- getLikelihoods(cD, nullData = TRUE, cl = cl)
The 'nullData = TRUE' option here is the bit that will treat separately those genes that show equivalent expression across timepoints, leaving those with high likelihood in the NDE group those that do not.
topCounts(cD, group = "NDE")
However, your model suggests that you believe there are two possible patterns to describe the data; one in which we have equivalent expression across all timepoints, and one in which "d2 and d5 should cluster together and differently from d0 and d60". To specifically model this, you'd need to modify the distributional assumptions in the densityFunction slot of the cD object. This is possible, but may be overkill - let me know if you'd like to see a discussion of this.