Hi!

I use the function dba.report to retrieve differentially bound sites (th = 1) I found the fold-changes tend to be very small and do not know how to compute them.

For example, at one site the mean for control is 1.6973 while the mean for treatment is 4.231, and the Fold is -0.001057009, p-value is 0.0051515283, FDR = 0.99.

How the fold-changes were computed? my code is below.

Thank you!

cond = "treatment"
obj = dba(sampleSheet=data0)
obj = dba.count(obj,minOverlap=2)
obj = dba.normalize(obj,normalize=DBA_NORM_RLE,library=DBA_LIBSIZE_PEAKREADS)
obj = dba.contrast(obj,reorderMeta=list(Condition=cond),minMembers=2)
obj = dba.analyze(obj)
dbp = dba.report(obj,th=1)

As you note, the fold change calculations are being performed by DESeq2, so you may want to look at the documentation for the DESeq2::lfcShrink function (and the References it lists) to see how it works.

My guess is that you are getting very low shrunken lfc values, and very high FDR values, due to having only two replicates which may not agree very well, leaving you with an under-powered experiment. You may want to check the read counts (by calling dba.report() with bCounts=TRUE) to see how consistent they are within sample groups.

I'm also wondering if you had a good reason to choose to normalize with library=DBA_LIBSIZE_PEAKS (RiP) as this tends to over-normalize, which also lessens apparent fold changes.