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"))
However, I noticed several strange things:
- 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.
- When I change the reference level to E, I do not get these significant genes.
- When testing C vs. D, I get different results depending on whether A or E are the refenece level.
- 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