I am testing DESeq2 1.10.1 version and find some issue with cleanContrast function.
I use a standard 2x2 A*B design with interaction terms.
DESeqDataSetFromMatrix(countData = round(datacounts), colData = designtable, design = ~ Gender*Site)
dds<- DESeq(dds, modelMatrixType="standard")
output <- results(dds, independentFiltering=T, contrast=c(0, 0, 1, 0))
When I extract the contrast results to get the site M vs site K when gender is Female with the contrast vector, I found that some of the all zero contrasts are not properly filtered. I debug the code using RStudio and find that in CleanContrast function , when a numeric contrast vector is used, contrastAllZero <- contrastAllZeroNumeric(object, contrast) is called. If I get the logic right, this function is used to detect genes whose counts for samples used in contrasts are all zero. However, after checking the source code of contrastAllZeroNumeric, it will do nothing for my contrast vector above.
### in contrastAllZeroNumeric function ###
if (all(contrast >= 0) | all(contrast <= 0)) {
return(rep(FALSE, nrow(object)))
}
I kind of feel this is a bug.
I should add, it works under the default DESeq() use of expanded model matrices, but not for standard model matrices.
Yeah, I think this function is written for expanded matrix only. For standard matrix, the filter works file when a contrast is specified in text as contrast=c("Term", "A", "B")
It makes biological sense to reset the log2fc and p-values for RNA-Seq data when all samples in the contrast have zero counts. Probably should extend this to complex models as well.
Thanks.