Hello All,
I'm working with a ChIP-seq dataset of RNA Pol II under two treatment conditions (NT and treated). I am interested in getting (and subsequently annotating) three peaksets, and these can be imagined as a Venn diagram of two circles: on the left is peaks that are DB in the NT, on the right is peaks that are DB in the treated, and in the middle is peaks that are occupied in both conditions but not DB. You can imagine that the peaks in that middle region would contain, among many other things, genes like Gapdh or Actb.
Acquiring the DB peaks was easy, as I just needed to do dba.report() with bGain or bLoss = TRUE. On the other hand, getting that middle peakset has been surprisingly difficult. If I run dba.report() with bNotDB = TRUE, I get a list of ALL peaks that have occupancy in at least one sample but don't meet the DB threshold. I managed to get a set of shared peaks by running dba.overlap() with mode = "DBA_OLAP_PEAKS", and the result$inAll peakset was what I was looking for. I could then subtract out the NotDB peakset from the inAll peakset to make that middle Venn diagram region.
However, DBA_OLAP_PEAKS only works with a maximum of four peaksets. Since I have 2 conditions with 2 replicates each, this works out today; however, if I run an experiment with more replicates, I need a new solution. Running dba.overlap() with mode = DBA_OLAP_RATE does not have a limit to number of samples, so the program is clearly able to perform the calculation. However, the output of that mode is simply a number, not a peakset that I can annotate.
How would you go about getting this shared peakset among all samples? I have been playing around with the consensus elements in the program, but it still doesn't make a ton of sense to me. Surely I'm overlooking something in DiffBind that would make this easy? I'm self taught and have been doing bioinformatics for less than a year, so it's quite likely I missed something simple.
Thank you!
