Hello,

I am using DiffBind to find differential binding sites for different ChIP-seq proteins. The experimental workflow consists to treat both wild type and transgenic mice to three different treatments: treatment 1, treatment 2 , treatment 3 (both T1+T2).  I wonder how it will be the best way to analyse the data in DiffBind. I saw that I can introduce a column in the csv table for the treatments and maybe another one for being wild type or transgenic, maybe also I can use this mask option described in the documentation, but I will really appreciate suggestions! 

Thanks!

S.

It depends on the question you want to ask. Are you most interested in differential reaction to a treatment by WT and transgenic mice? Are you looking for how the reactions differ between treatment, independently of whether they are WT or transgenic? 

Let's assume your DBA object is called myDBA, and you use DBA_CONDITON to indicate WT or Transgenic, and DBA_TREATMENT to be T1, T2, or T1+T2.

To identify deferentially bound sites between Treatment 1 and Treatment 2 independent of mouse type, your could use:

> myDBA = dba.contrast(myDBA, group1=myDBA$masks$T1, group2=myDBA$masks$T2, block=DBA_CONDITION, name1="T1", name2="T2")

To identify sites that are deferentially bound in WT vs Transgenic mice for Treatment 1:

> myDBA = dba.contrast(myDBA, group1=myDBA$masks$T1 & myDBA$masks$WT, group2=myDBA$masks$T1 & myDBA$masks$Transgenic, name1="WT:T1", name2="Transgenic:T1")

Cheers-

Rory

Thanks a lot!! 

Sergio

Ok, I understand now much better, I guess that I can not answer everything at the same time. If I want to asses the differential binding sites of each treatment I need to perform different contrast right? 

If I want only to asses the differential binding sites only of T1 compared to untreated samples, could I proceed this way?:

> myDBA = dba.contrast(myDBA, group1=myDBA$masks$Untreated, group2=myDBA$masks$T1, block=DBA_CONDITION, name1="Untreated", name2="T1")

Thanks a lot Rory

If I want to asses the overlapping between:

> myDBA = dba.contrast(myDBA, group1=myDBA$masks$T1, group2=myDBA$masks$T2, block=DBA_CONDITION, name1="T1", name2="T2")

> myDBA = dba.contrast(myDBA, group1=myDBA$masks$T1 & myDBA$masks$WT, group2=myDBA$masks$T1 & myDBA$masks$Transgenic, name1="WT:T1", name2="Transgenic:T1")

Should I proceed directly with :

dba.overlap(myDBA, group1=myDBA$masks$T1, group2=myDBA$masks$T2, mode=DBA_OLAP_RATE)

dba.overlap(myDBA, group1=myDBA$masks$T1 & myDBA$masks$WT, group2=myDBA$masks$T1 & myDBA$masks$Transgenic, name1="WT:T1", name2="Transgenic:T1", mode=DBA_OLAP_RATE)

Sorry , maybe it is really naive, but I think that I should maybe generate two different objects? I am a bit lost. 
 

Hi Sergio-

I think maybe there is some confusion between two styles of analysis.

dba.overlap() and dba.plotVenn() are used for an occupancy analysis, which identifies how called peaks relate (which ones are unique to one sample and which ones overlap). dba.overlap() does not deal with sample groups -- if you have replicates, you should create  consensus peaksets (using dba.peakset()) first and then see how they overlap.

dba.contrast() and dba.analyze() are used for an affinity analysis, which looks quantitatively at read counts and identifies regions where the binding strength changes consistently between replicated sample groups using a differential expression package (edgeR or DESeq2).

You can set up all the contrasts you want by calling dba.contrast() for each contrast (or using its ability to generate a bunch of contrasts at once), then a single call to dba.analyze().  You can do this all within a single DBA object. For each contrast, calling dba.report(contrast=n) will return the sites identified as being significantly differentially bound in contrast n, sorted by FDR value (lowest first).

Cheers-

Rory

Ok, sorry , I see the issue.

I thought that I could compare the the different DBS called by different affinity analysis doing a Venn Graph. In theory I could do one to see if using different affinity analysis can share the same overlapping regions or not , right? But I guess I will need to create a peak set of DBS identified for each analysis.  I think could be nice to see synergies among different treatments 

Hi Sergio-

The easiest way to do this is using dba.report() to generate a new DBA object containing the DB (gain and loss) peaksets.  The man page for dba.report() (obtained by typing  ?dba.report) has an example of how to do what you want. It is the last example on the page, headed #Retrieve report-based DBA object, as follows:

#Retrieve report-based DBA object
data(tamoxifen_counts)
tamoxifen <- dba.contrast(tamoxifen, categories=DBA_CONDITION, block=tamoxifen$masks$MCF7)
tamoxifen <- dba.analyze(tamoxifen,bCorPlot=FALSE)
tamoxifen.DB <- dba.report(tamoxifen,method=c(DBA_EDGER,DBA_EDGER_BLOCK),
                          bDB=TRUE, bGain=TRUE, bLoss=TRUE)
dba.plotVenn(tamoxifen.DB,1:4,label1="Single Factor GAIN",label2="Single Factor LOSS",
                             label3="Blocking Factor GAIN",label4="Blocking Factor LOSS")

Cheers-

Rory