Greetings,

I am using DiffBind for a differential chromatin accessibility analysis.

During manipulation of the output DARs, I noticed an inconsistency regarding the number of the exported regions:

# using dba.report filtering parameters
ATAC.dup.DESEQ2.DB.fdr <- dba.report(ATAC.dup.DESEQ2.analyze,method=DBA_DESEQ2,th=0.1,bUsePval=FALSE,fold=1.5)
df <-as(ATAC.dup.DESEQ2.DB.fdr, "data.frame")
nrow(df)
[1] 355

while when I perform the filtering manually:

# using custom filtering parameters
ATAC.dup.DESEQ2.DB.fdr <- dba.report(ATAC.dup.DESEQ2.analyze,method=DBA_DESEQ2,th=1,bUsePval=FALSE,fold=0)
df <-as(ATAC.dup.DESEQ2.DB.fdr, "data.frame")
nrow(df[which(df$FDR<=0.1 & abs(df$Fold)>=1.5),])
[1] 368

So, when I use the dba.report filter parameters I actually lose 13 differentially accessible regions:

seqnames     start       end width strand Conc Conc_UV2h Conc_NOUV  Fold        p.value      FDR
chr3    130679103    130679763    661    *    6.51    7.08    5.58    1.5    1.77e-06    0.000683
chr2    182410605    182411350    746    *    6.35    6.91    5.41    1.5    2.15e-06    0.000722
chr12    69921416    69922192    777    *    6.2    6.77    5.27    1.5    2.78e-05    0.00151
chr18    57213269    57213871    603    *    5.73    6.29    4.79    1.5    0.000103    0.00236
chr4    28371089    28371630    542    *    5.24    5.8    4.31    1.5    0.00106    0.0074
chr1    101232858    101233451    594    *    4.95    5.52    4.02    1.5    0.00181    0.0103
chr2    67495854    67496408    555    *    4.81    5.38    3.88    1.5    0.00214    0.0114
chr12    29061635    29062631    997    *    4.73    5.29    3.8    1.5    0.00258    0.0129
chr7    152038356    152038917    562    *    4.73    5.29    3.8    1.5    0.00272    0.0133
chr6    131494779    131495406    628    *    4.82    5.38    3.88    1.5    0.00393    0.0169
chr17    30788814    30789641    828    *    4.49    5.05    3.56    1.5    0.0049    0.0197
chr12    28152380    28153031    652    *    4.57    5.14    3.64    1.5    0.00562    0.0216
chr1    165570994    165571154    161    *    4.5    5.06    3.56    1.5    0.00674    0.0244

Also note that dba.report help page says: "..less/greater than or equal to this value will be included to the report.." for both th and fold parameters, so using "<= "and ">= " in the custom filtering command is valid. 

So I think that there is a probable bug here, or I am missing something.

Cheers,
Dimitris

Hello Dimitris-

The issue you raise is related to the precision of the calculations. The returned report has reduced precision for the Fold values (only one place to the right of the decimal point), so these 13 sites actually have log Fold values slightly less than 1.5 and are being rounded to 1.5 in the report. Within dba.report(), the filtering is done on full-precision values, so these sites don't meet the Fold cutoff criterion.

So you can either use the reporting function to do the filtering at full precision, or filter the report yourself at lower precision to "rescue" these intervals. I may consider adding the desired precision as a parameter to dba.report() in future for finer grain control over this.

Hope this helps-

Rory

Hello Rory,

ok, your answer is very clear.

Thank you,

Dimitris