Question

Trying to use equivalence test from a diffbind dataset

2

Entering edit mode

6 months ago

jandrade ▴ 20

Hello,

I have an ATAC-seq dataset that I've analyzed using Diffbind. So far so good, but one of the things I'm interested in is whether some peaks that change from time point B to C return to a state similar to an earlier time point A. I'm aware that DESeq2 itself has the option to run a test of equivalence using

results(dds, lfcThreshold=lfcThreshold, altHypothesis="lessAbs")

I think there also might be an option using

lfcShrink(dds, coef=coef_name, type="apeglm")

The issue that I'm having though is pulling out the DESeq2 object to then use with DESeq2 directly. My DBA object is already analyzed and has the right number of contrasts.

> dba.show(data, bContrasts = T)
      Factor Group Samples Group2 Samples2 DB.DESeq2
1     Factor  bio1      17   bio2       14     68356
2  Condition 26hpf       5  42hpf        6     56535
3  Condition 26hpf       5  72hpf        5     73584
4  Condition 26hpf       5   7dpf        5     83941
5  Condition 26hpf       5   1hpb        6     84790
6  Condition 26hpf       5   3dpb        4     79843
7  Condition 42hpf       6  72hpf        5     48889
8  Condition 42hpf       6   7dpf        5     77501
9  Condition 42hpf       6   1hpb        6     78192
10 Condition 42hpf       6   3dpb        4     68424
11 Condition 72hpf       5   7dpf        5     61107
12 Condition 72hpf       5   1hpb        6     59856
13 Condition 72hpf       5   3dpb        4     53611
14 Condition  7dpf       5   1hpb        6      9281
15 Condition  7dpf       5   3dpb        4     45207
16 Condition  3dpb       4   1hpb        6     48542

I tried to extract the deseqdataset object using these two methods, which both give me the same thing

dds <- dba.analyze(data, bRetrieveAnalysis = T)
dds <- data$DESeq2$DEdata

But I get the wrong number of contrasts.

> resultsNames(dds)
[1] "Intercept"                "Factor_bio2_vs_bio1"     
[3] "Condition_42hpf_vs_26hpf" "Condition_72hpf_vs_26hpf"
[5] "Condition_7dpf_vs_26hpf"  "Condition_1hpb_vs_26hpf" 
[7] "Condition_3dpb_vs_26hpf"

Is this a bug? How can I make sure all of my contrasts get included? That's problem 1

Problem 2 comes from the actual analysis. I'm just not sure if I'm doing this correctly. I tried to plug just one of my contrasts into DESeq2::results() and got some NAs in the pvalue column.

> res_eq <- results(dds, altHypothesis="lessAbs", lfcThreshold = log2(1.25))
> res_eq[res_eq$pvalue<0.05,]
Error: logical subscript contains NAs
> res_eq[which(is.na(res_eq$pvalue)),]
log2 fold change (MLE): Condition 3dpb vs 26hpf 
Wald test p-value: Condition 3dpb vs 26hpf 
DataFrame with 1 row and 6 columns
       baseMean log2FoldChange     lfcSE      stat    pvalue
      <numeric>      <numeric> <numeric> <numeric> <numeric>
29674   23.4258       -2.21303  0.759665         0        NA
           padj
      <numeric>
29674        NA

I have no idea why these show up and it makes me think I'm doing something wrong.

Any insight would be appreciated, whether that's on extracting the deseq data object correctly, running the equivalence test correctly, or just on whether there's a better way to answer my question. Thanks!

Diffbind R DESeq2 ChIP-seq ATAC-seq • 489 views

ADD COMMENT • link updated 6 months ago by ATpoint 86k • written 6 months ago by jandrade ▴ 20

score 2 · Answer 1 · 2024-06-13

Is this a bug? How can I make sure all of my contrasts get included? That's problem 1

apeglm is cumbersome, you might need to relevel several times to test all possible contrasts: type='apeglm' shrinkage only for use with 'coef'

I have no idea why these show up and it makes me think I'm doing something wrong.

See FAQs: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#why-are-some-p-values-set-to-na

Also, be sure to read the DESeq2 manual on contrasts.

res_eq <- results(dds, altHypothesis="lessAbs", lfcThreshold = log2(1.25))

By experience, the lessAbs needs quite some power, and using 1.25 might be was to stringent to get good results. You might be forced to use something like 1.5 to get a good number of genes. Otherwise you really just enrich for genes with large counts that have a lot of power and are not DE.