Discrepancy in DE results from base DESeq2 and tidybulk DESeq2 method

[MSaadfarooq]

Hi, I used following code for DESeq2 workflow but found some differences in results

gse
# A SummarizedExperiment-tibble abstraction: 3,532,606 × 30
# Features=36047 | Samples=98

dds <- DESeqDataSet(gse, design = ~ batch + sex + status) 
dds <- DESeq(dds)
res <- results(dds)
summary(res)
out of 35983 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 188, 0.52%
LFC < 0 (down)     : 62, 0.17%
outliers [1]       : 99, 0.28%
low counts [2]     : 15312, 43%
(mean count < 7)

I also performed same SE object with same variables using tidybulk

da <- gse |> test_differential_abundance(~batch + sex + status, method = "DESeq2")
da_res <- da %>% pivot_transcript()

# to compare results
all.equal(rownames(res), da_res$.feature)
[1] TRUE
sum(da_res$padj < 0.1, na.rm=TRUE)
[1] 321
table(base_sig = res$padj < .1, tidy_sig = da_res$padj < .1)
        tidy_sig
base_sig FALSE  TRUE
   FALSE 20119   255
   TRUE    224    26

Please guide where it went wrong?

[stemangiola]

@mikelove FYI

[mikelove]

See recent support site post on same question... scaling needs to be set also.

[stemangiola]

@MSaadfarooq do you mind crossreferencing the bioconductor issue here?

[MSaadfarooq]

Please go ahead, will check that, Thank you!

[stemangiola]

@zhanchen can you please take care of this? It will be a easy fix in the DESeq2 stream , for Summarized experiment methods as well as tibble methods.

[zhanchen]

Hey @MSaadfarooq, @stemangiola @mikelove
I had a look at the code and ran some tests, and the differences between above results come from arg contrast is not set in your code res <- results(dds) for the function results. However, tidybulk will automatically set contrast in its code https://github.com/stemangiola/tidybulk/blob/dcbf6598fbf761b4bba068c7aee07cf9ae9f4eec/R/functions_SE.R#L1286C6-L1290C7

Arg contrast from function results has no default value but can accept missing arg. Also, function test_differential_abundance listed arg scaling_method = "TMM", but seems did nothing related to scaling.

[mikelove]

scaling_method should be RLE to be similar to DESeq2's median ratio method.

[stemangiola]

We should definetly use the DESeq2 normalisation, (passed through the scaling_method argument) we need 100% concordance. I thought tidybulk was doing that already. But if DESeq2 has only one scaling method, then we just use that.

Accordingly, we should improve the documentation, relative to this topic for DESeq2