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?

I think it is because the default normalization is not DESeq2.

Try setting the scaling_method within test_differential_abundance(). You can further explore this with:

scale_abundance(method="RLE") # DESeq2 scaling

There may still be minor differences, as this is the edgeR implementation of DESeq2 scaling.

Also you should turn off independentFiltering in DESeq2 when calling results() if you want to compare 1-to-1. You can filter by abundance with tidybulk at the top.