I am a little puzzled by the way of generating a dba report from DiffBind matrix. When I run dba.report(mat, contrast=1), this should run with default fold and p-values, which should be 1 (fold) and 0.05 (p-value). Then, if I "relax" the fold parameter [dba.report(mat, contrast=1, fold=0.5)], theoretically, I should get the same features as with fold=1 but with additional ranges in the results. What is more puzzling is that I "additionally" subset the resulting GRanges using the Fold column to select only the positive fold change ranges and exclude the negative ones [for example, using mat.DB1$Fold>1,]. Again, in theory, I should be able to get the same GRanges in both approaches, but I am not.

Is the default value in dba.report(mat, contrast=1) for fold cut off not 1? Or is dba.report(mat, contrast=1, fold=0.5) doing something different from what I am expecting?

I came to identify this discrepancy because I noticed one locus has a Fold of 1.7311, it only appears in my GRanges when I use fold=0.5 and not when I use default value, despite the p-value being so small (10^-16), and should have been met the cut off with both approaches.

Thanks for your time.

dba.report(mat, contrast=1)
# GRanges object with 21960 ranges and 6 metadata columns
# but with specifying fold=0.5
dba.report(mat, contrast=1, fold=0.5)
# GRanges object with 1261 ranges and 6 metadata columns:


sessionInfo( )

When you specify a fold threshold to dba.report() the null hypothesis is changed and all of the statistics are re-computed, so each interval has a different FDR value.. This is documented for the fold parameter on the help page for dba.report().

You can see how the FDR values change by comparing the full reports:

report1 <- dba.report(mat, contrast=1, fold=0.0)
report2 <- dba.report(mat, contrast=1, fold=0.5)