Question: DESeq2: name and contrast arguments return partially different results under specific circumstances -- bug?
0
4 months ago by
Vincent de Bakker0 wrote:

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

deseq2 • 118 views
modified 4 months ago by Michael Love26k • written 4 months ago by Vincent de Bakker0
Answer: DESeq2: name and contrast arguments return partially different results under spe
1
4 months ago by
Michael Love26k
United States
Michael Love26k wrote:

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.

Ok thanks for the prompt answer.

It would indeed be great if you could mention this clearly in the vignette! This may in these very specific cases seriously alter what researchers will report as up- or down-regulated. The difference in reported LFCs that got me looking into this was huge (so I knew something was going on). In my case, the few different LFCs I obtained using name were just really not right. I searched the web long and thoroughly to find more information about it, but found none.

I'll just use the contrast argument from now on by default whenever I can.

The LFCs can be large without moderation, which is the default with DESeq(), we recommend if you are looking at the LFC to apply shrinkage via lfcShrink, which would shrink unreliably LFC as well.