Hi,

In the DESeq2 vignette, it is stated specifically that the name and contrast argument in results() can be used to build equivalent results tables (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#differential-expression-analysis). I found this in general to be true, with one exception. I think this might be a bug.

Explanation

The results are only different when:

1. We have more than two conditions in our DESeqDataSet, and
2. the gene has exclusively zero counts in both conditions to be contrasted, but not in at least one sample of another condition.

If the gene had only zeroes in all samples, we would get NA values. We do not get that here, because the nonzero counts in the other condition allow still for the estimation of the dispersion parameter.

Below I provide a minimal example. Although the difference in output of that example seems marginal, I found it can be very big with real data. I have seen differences in LFC of -22.5 (using name) and 0 (using contrast), for the exact same data and contrasts. The corresponding p-values can then have huge differences too.

Questions

1. Is this indeed a bug, or do I misinterpret something?
2. Are the results of DESeq2 obtained using the contrast argument therefore more trustworthy?
3. In case this is not a bug, could anyone explain to me where the difference comes from? Are these two calls not supposed to perform exactly the same computations and tests?
4. In case yes to question 2: is there a way to extract interaction term results using the contrast argument?

Minimal, reproducible working example

library(DESeq2)
# reproducible minimal data set 
set.seed(1992)
dds <- makeExampleDESeqDataSet(m = 9)
colData(dds)$condition <- factor(rep(LETTERS[1:3], each = 3)) # 3 conditions
colData(dds)

So we have a simulated data set with 1000 genes, 9 samples, 3 conditions and 3 samples per condition:

# output:
DataFrame with 9 rows and 1 column
        condition
         <factor>
sample1         A
sample2         A
sample3         A
sample4         B
sample5         B
sample6         B
sample7         C
sample8         C
sample9         C

We run DESeq and build the results table contrasting conditions A and B in the next chunk. We can do this by using either the name or the contrast argument. Later on we would be interested in the other contrasts with condition C as well, but for purposes of illustration I leave that out in this example.

dds <- DESeq(dds)
# call results() w/ name and contrast arguments
res.cont <- results(dds, contrast = c("condition", "B", "A"), lfcThreshold = 1)
res.name <- results(dds, name = "condition_B_vs_A", lfcThreshold = 1)
# at first glance same results
res.cont; res.name

Indeed, the results returned appear to be the same:

# ouput using contrast
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 1000 rows and 6 columns
                 baseMean     log2FoldChange             lfcSE               stat            pvalue      padj
                <numeric>          <numeric>         <numeric>          <numeric>         <numeric> <numeric>
gene1    24.9730902462296  0.513386216808375 0.694661298745144                  0                 1         1
gene2    44.8672687070342 -0.141742666322917 0.491130886311999                  0                 1         1
gene3    13.1340089761356 -0.524733081118937 0.840581037100424                  0                 1         1
gene4    4.90353109443876  -2.52272107359962  1.34233233414978  -1.13438455951677 0.256633272980206         1
gene5    22.5591506152625 -0.518806678226352 0.702241365250218                  0                 1         1
...                   ...                ...               ...                ...               ...       ...
gene996  2.98439046433484  0.978545031361552  1.84169185147413                  0                 1         1
gene997  7.15961550810323 -0.233795493067082 0.949957082688798                  0                 1         1
gene998  52.7301561370165  0.599617677250584 0.496828957431585                  0                 1         1
gene999  30.8309983074886 -0.246064160021411 0.592298347944716                  0                 1         1
gene1000 12.3788663640862  -1.32340176110781 0.887650657934792 -0.364334502787657 0.715608238264671         1

Above the output using contrast, below for name:

# output using name
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 1000 rows and 6 columns
                 baseMean     log2FoldChange             lfcSE               stat            pvalue      padj
                <numeric>          <numeric>         <numeric>          <numeric>         <numeric> <numeric>
gene1    24.9730902462296  0.513386216808375 0.694661298745144                  0                 1         1
gene2    44.8672687070342 -0.141742666322917 0.491130886311999                  0                 1         1
gene3    13.1340089761356 -0.524733081118937 0.840581037100424                  0                 1         1
gene4    4.90353109443876  -2.52272107359962  1.34233233414978  -1.13438455951677 0.256633272980206         1
gene5    22.5591506152625 -0.518806678226352 0.702241365250218                  0                 1         1
...                   ...                ...               ...                ...               ...       ...
gene996  2.98439046433484  0.978545031361552  1.84169185147413                  0                 1         1
gene997  7.15961550810323 -0.233795493067082 0.949957082688798                  0                 1         1
gene998  52.7301561370165  0.599617677250584 0.496828957431585                  0                 1         1
gene999  30.8309983074886 -0.246064160021411 0.592298347944716                  0                 1         1
gene1000 12.3788663640862  -1.32340176110781 0.887650657934792 -0.364334502787657 0.715608238264671         1

However, there are a few genes for which the LFC values are not the same:

# find for which transcripts output is not equal
ind <- which(res.cont$log2FoldChange != res.name$log2FoldChange)
res.cont[ind, ]; res.name[ind, ]

For contrast:

# output contrast
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 3 rows and 6 columns
                 baseMean log2FoldChange            lfcSE      stat    pvalue      padj
                <numeric>      <numeric>        <numeric> <numeric> <numeric> <numeric>
gene95   1.56000455665874              0 2.94271443607031         0         1         1
gene198 0.324896579208153              0  4.0804558717943         0         1         1
gene227 0.102925664787549              0 4.08045587181216         0         1         1

...and for name:

# output name
log2 fold change (MLE): condition B vs A 
Wald test p-value: condition B vs A 
DataFrame with 3 rows and 6 columns
                 baseMean      log2FoldChange            lfcSE      stat    pvalue      padj
                <numeric>           <numeric>        <numeric> <numeric> <numeric> <numeric>
gene95   1.56000455665874 -0.0473928419930928 2.94271443607031         0         1         1
gene198 0.324896579208153 -0.0474023370178275  4.0804558717943         0         1         1
gene227 0.102925664787549 -0.0473898765006612 4.08045587181216         0         1         1

When we check the (normalized) counts of these specific genes, we see the pattern I described above:

# check normalized counts of those transcripts
counts(dds, normalized = TRUE)[ind, ]

The counts for all samples we compare here are zero, but not so for all remaining samples in condition C:

# output
        sample1 sample2 sample3 sample4 sample5 sample6   sample7 sample8  sample9
gene95        0       0       0       0       0       0 11.115972       0 2.924069
gene198       0       0       0       0       0       0  0.000000       0 2.924069
gene227       0       0       0       0       0       0  0.926331       0 0.000000

The answer we look for is an LFC of zero, since there is no difference between the samples we compare. This result is given by a call to results() using contrast, but not using name.

Session info

> sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq2_1.22.2               SummarizedExperiment_1.12.0 DelayedArray_0.8.0          BiocParallel_1.16.6         matrixStats_0.54.0          Biobase_2.42.0              GenomicRanges_1.34.0       
 [8] GenomeInfoDb_1.18.2         IRanges_2.16.0              S4Vectors_0.20.1            BiocGenerics_0.28.0        

loaded via a namespace (and not attached):
 [1] locfit_1.5-9.1         Rcpp_1.0.1             lattice_0.20-38        zeallot_0.1.0          digest_0.6.20          backports_1.1.4        acepack_1.4.1          RSQLite_2.1.1         
 [9] ggplot2_3.2.0          pillar_1.4.2           zlibbioc_1.28.0        rlang_0.4.0            lazyeval_0.2.2         rstudioapi_0.10        data.table_1.12.2      annotate_1.60.1       
[17] blob_1.2.0             rpart_4.1-13           Matrix_1.2-15          checkmate_1.9.4        splines_3.5.2          geneplotter_1.60.0     stringr_1.4.0          foreign_0.8-71        
[25] htmlwidgets_1.3        RCurl_1.95-4.12        bit_1.1-14             munsell_0.5.0          compiler_3.5.2         xfun_0.8               pkgconfig_2.0.2        base64enc_0.1-3       
[33] htmltools_0.3.6        nnet_7.3-12            tibble_2.1.3           gridExtra_2.3          htmlTable_1.13.1       GenomeInfoDbData_1.2.0 Hmisc_4.2-0            XML_3.98-1.20         
[41] crayon_1.3.4           bitops_1.0-6           grid_3.5.2             xtable_1.8-4           gtable_0.3.0           DBI_1.0.0              magrittr_1.5           scales_1.0.0          
[49] stringi_1.4.3          XVector_0.22.0         genefilter_1.64.0      latticeExtra_0.6-28    vctrs_0.2.0            Formula_1.2-3          RColorBrewer_1.1-2     tools_3.5.2           
[57] bit64_0.9-7            survival_2.43-3        yaml_2.2.0             AnnotationDbi_1.44.0   colorspace_1.4-1       cluster_2.0.7-1        memoise_1.1.0          knitr_1.23            

This is known behavior (described in NEWS sections for 1.21 and 1.15), but I can certainly add more description to the vignette, so it's better understood.

Using contrast, we zero out LFCs when the two groups specified by the contrast have all zero counts. It's non-trivial to do this by name, and so it's only a feature for the contrast.