Hi Michael,

I have couple of questions

1. The number of differentially expressed genes shot almost 6 fold (from ~1000 to ~6000) between deseq and deseq2. Is this a common trend and if so, why?

2. For gene, the raw counts are as follows (control - 0,10,670; treatment - 12986,9118, 6409). But the results looks like this. NAs are generated because of one zero value in the control? I thought pseudo counts are added if the counts are zero.

 baseMean        log2FoldChange     lfcSE           stat                pvalue       padj
3331.195751      4.6590639             1.2223685   3.8115050          NA         NA

when I converted the ZEROs to ONEs pvalues are as follows...and the number of differentially expressed genes also reduced.

4076.763438     4.29728935           0.9134027   4.70470420     2.542343e-06  6.884325e-05

hi Prasad,

1. Yes this is expected. See the manuscript (linked from the first page of the vignette) for a discussion of the difference in dispersion estimation, the use of independent filtering to increase power, and the increased sensitivity in the results from DESeq to DESeq2. Also in our manuscript we discuss the situation when there are many significant differences, and how to find genes with significantly large effect sizes (as opposed to nonzero effect sizes). You can reduce the size of the list you are interested in by either lowering the alpha or using the lfcThreshold argument of results().

2. The NA are discussed in the vignette (see the FAQ which has a link to the section). With 3 samples per group, the filtering on Cook's distance might not be ideal for your dataset. You can turn it off and investigate these genes manually.

We do not recommend altering the counts.

Thank you Michael. I set minReplicatesForReplace=Inf and reran the DESeq but this didn't change the number of differentially expressed genes (Same with the default number for minReplicatesForReplace).

Also I set cooksCutoff to FALSE and increased the number of differential genes. 

Here is the sample code

dds <- DESeqDataSetFromMatrix(countData = countsMatrix, colData = colData, design = ~ type);
dds <- DESeq(dds);

Gres <- results(dds, contrast=c("type","ABCD_DIF","ABCD_UND"), cooksCutoff = FALSE);