Clarification on contrast levels in DESeq2

Question

[Menachem Sklarz]

Hi,

I was running DESeq2 with the simplest possible design: one factor (called group) with 5 levels: A-E, where A is the reference level.

I run the model thus:

dds <- DESeqDataSetFromMatrix(countData = counts_data,
                              colData = col_data,
                              rowData = row_data,
                              design = ~ group)
dds <- DESeq(dds)

To test pairwise contrasts, I do the following:

b<-results(dds,contrast = c("group","A","B"))

etc.

However, I noticed several strange things:

1. For the A-B comparison, I get significant genes when all samples in groups A and B have zero counts! It looks as if the comparison is against the Intercept, rather than level A, even though the results object says the data is for A vs. B comparison.
2. When I change the reference level to E, I do not get these significant genes.
3. When testing C vs. D, I get different results depending on whether A or E are the refenece level.
4. I get a significant enrichment (adjusted p-value of 1.06E-04) when the reads in group A are 0 (for two samples) and the reads in group B are c(1,2,2). Does this make sense? Aren’t extremely low-count genes supposed to be filtered out in the “independent filtering” stage?

Note: I solved issues 1-3 by setting betaPrior to TRUE and modelMatrixType to "expanded".

Thanks a lot! Menachem

[Kevin Blighe]

An obvious comment: why are you even using data for samples that are all zero?