I'm trying to do "combine variables" approach in a two-factor experimental design. To create combinations of variables I follow the example from here http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions

Strangely, the list of differentially expressed genes is different based on how I setup colData:

strain <- c('WT', 'WT', 'KO', 'KO', 'KO', 'WT','WT','WT','WT' ,'KO', 'KO','KO')
condition <- c('A','A','A','A','A','B','B','B','B','B','B')
colData <- data.frame(sample = colnames(counts),
                      group = factor(paste0(strain, "_", condition),
                              levels=c('WT_A', 'WT_B', 'KO_A', 'KO_B')),
                      stringsAsFactors = FALSE
) 
#same thing without explicitly listing factor levels
colData <- data.frame(sample = colnames(counts),
                      group = factor(paste0(strain, "_", condition)                                  
                      stringsAsFactors = FALSE
) 

#these two colData objects produce different results when I run
dds <- DESeqDataSetFromMatrix(countData = counts[, colData$sample],
                              colData = colData,
                              design =~ group)
dds <- dds[rowMeans(counts(dds)) >= 10, ]
dds <- DESeq(dds)
result <- results(dds, contrast = c('group', 'WT_B', 'WT_A'), 
                       cooksCutoff = FALSE, 
                       independentFiltering = FALSE)

Interestingly, the first approach gives the same result as a classic linear model: strain + condition + strain:condition. Is there some subtle issue with how DESeq expects the factors to be supplied in the design?

The difference in DE gene lists is also quite peculiar. It is not like all genes are scrambled, rather few of them go up and down the rank. Here is the top 10 genes and only one of them is misplaced in each set.

Irx2 Zic2 Ltbp1 Etv4 Akr1b7 Mgam Efhd1 Ugt1a6a Sgms2

Irx2 C1s1 Ltbp1 Etv4 Akr1b7 Mgam Efhd1 Ugt1a6a Sgms2

Can you take a look at a plot of the p-values (you don't need to post)? Is it the case that the p-values are the same except for a few genes? Or are the p-values different for most or all genes?