Hi Rory,

I was wondering whether Diffbind could call differential peaks from bed file? The reason I ask this is that if I would like to all differential peaks from two sets of diffPeaks called by diffBind.

For example, I have group1, group2, group3, group4.

Then I already used diffbind to call differential Peaks. Diffpeak1=diffbind(group1,group2), Diffpeak2=diffbind(group3,group4).

Then Is there a way that I could get diffbind(Diffpeak1,Diffpeak2).

Thank you so much.

Y

Hello Yijing-

I'm going to assume that your DBA_CONDITION values are cond1 and cond2, and your DBA_TREATMENT values are treat1 and treat2. You are first contrasting cond1:treat1 and cond1:treat2 (getting Diffpeak1) and cond2:treat1 and cond2:treat2 (getting Diffpeak2). Now you want to contrast cond1 with cond2, looking only at the peaks in the union of Diffpeak1 and Diffpeak2. Is that a resonable way of re-phrasing the question?

There area number of ways to do this in DiffBind. Here is one that should work::

> analysis <- dba.contrast(DBA,group1=DBA$masks$cond1&DBA$masks$treat1, name1="cond1:treat1",
                          group2=DBA$masks$cond1&DBA$masks$treat2, name2="cond1:treat2")

> analysis <- dba.contrast(analysis,group1=DBA$masks$cond2&DBA$masks$treat1, name1="cond2:treat1",
                          group2=DBA$masks$cond2&DBA$masks$treat2, name2="cond2:treat2")
> analysis <- dba.analyze(analysis)

​​Now you can get the differentially bound sites in a new DBA object using the report-based object feature of dba.report():

> DB <- dba.report(analysis,bDB=TRUE)
> DiffPeak1 <- dba.peakset(DB, peaks=1, bRetrieve=TRUE)
> DiffPeak2 <- dba.peakset(DB, peaks=2, bRetrieve=TRUE)
> overlaps <- dba.overlap(DB, 1:2)
> uniqueCond1 <- overlaps$onlyA
> uniqueCond2 <- overlaps$onlyB
> overlapping <- overlaps$inAll

At this point, overlapping would contain the peaks that were differential in both conditions (driven by the treatment). The two unique sets would contain peaks that were changed by the treatment, but only in one or the other of the conditions.

If you want to do another differential analysis, you can generate the a peakset with any peaks in either of these sets using the consensus feature of dba.peakset():

> DB <- dba.peakset(DB, consensus=TRUE, minOverlap=1)
> DBpeaks <- dba.peakset(DB, peaks=3, bRetrieve=TRUE)

Now you can use these peaks to re-count the original data and run another analysis:

> newDBA <- dba.count(DBA, peaks=DBpeaks, minOverlap=1)

I'm not sure exactly what groups you would compare here -- perhaps all the samples for each of the conditions? In that case your would probably want to using a blockg design to model the treatments:

> analysis <= dba.contrast(newDBA, group1=newDBA$masks$cond1, name1="cond1",
                           group2=newDBA$masks$cond2, name2="cond2"
                           block=DBA_TREATMENT)
> analysis <- dba.analyze(analysis)

The differentially bound sites from this would be ones that were consistently changed between by two conditions.

You could also run an analysis in one go, from the beginning:

> analysis <= dba.contrast(DBA, group1=DBA$masks$cond1, name1="cond1",
                           group2=DBA$masks$cond2, name2="cond2"
                           block=DBA_TREATMENT)

> analysis <- dba.analyze(analysis)

The resulting differentially bound sites would be ones that were consistently different between the two conditions (within treatment groups). 

Hope this helps -- you are trying to do something fairly complicated, and there may be better ways to do this using a more sophisticated design matrix than DiffBind can handle.

Cheers-

Rory

Hi Rory,

Thank you so much for your detail explanation.

The way you display my question is much better. I really appreciated your help. The strategies you provided here are very helpful. Thank you again.

Yijing