Question: DESeq2: name and contrast arguments return partially different results under specific circumstances -- bug?
0
gravatar for Vincent de Bakker
4 weeks 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 • 87 views
ADD COMMENTlink modified 4 weeks ago by Michael Love24k • written 4 weeks ago by Vincent de Bakker0
Answer: DESeq2: name and contrast arguments return partially different results under spe
1
gravatar for Michael Love
4 weeks ago by
Michael Love24k
United States
Michael Love24k 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.

ADD COMMENTlink written 4 weeks ago by Michael Love24k

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.

ADD REPLYlink written 4 weeks ago by Vincent de Bakker0

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.

ADD REPLYlink written 4 weeks ago by Michael Love24k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 288 users visited in the last hour