consensus analysis and dba.count function
4
0
Entering edit mode
pabla86 • 0
@pabla86-12186
Last seen 7.2 years ago

Dear author,

I am using your package Diffbind to perform chip-seq analysis of my data. In my analysis, I used three different peak callers for each sample, and I used Diffbind to retrieve the consensus. Now, I would like to center the consensus peaks in order to perform a motif analysis. I read (https://www.biostars.org/p/87925/) that  you say that Diffbind allows it by using the dba.count function and using the "summits" argument.

I tried to perform this analysis, here the the script:

sample1 = dba.peakset(sample1, minOverlap=3,sampID="sample1_consensus")

sample1<-dba.count(sample1,peaks=sample1$masks$Consensus,summits=200)


but the following error appears:
Error in pv$peaks[[i]] : index out of bound

I don't know how to solve this error. I would also like t ask you how I can retrieve the genomic coordinates of the centered peaks after dba.count function.
Thank you in advance for your help.
Best regards,

Paola

diffbind consensus • 1.7k views
ADD COMMENT
0
Entering edit mode
Rory Stark ★ 5.1k
@rory-stark-5741
Last seen 8 days ago
Cambridge, UK

Hi Paola-

What version of DiffBind are you using? 

There is a known bug that prevented recent version of DiffBind from computing summits when specifying a consensus peakset. Thus bug has been resolved in the more recent version, 2.2.7. If you can update to that, it should take care of the problem.

Another possible workaround in your case, which would also simplify the code you sent, would be to allow dba.count() to compute the consensus peakset implicitly:

> sample1 <- dba.peakset(sample1, minOverlap=3,sampID="sample1_consensus")
> sample1 <- dba.count(sample1,minOverlap=3,summits=200)

Hope this helps-

Rory

ADD COMMENT
0
Entering edit mode
pabla86 • 0
@pabla86-12186
Last seen 7.2 years ago

Dear Rory,

I 've tried both options, and they didn't work. In  particular, after updating DiffBind, and performing the analysis as you suggested, I've got the following error:

sample1 <-dba.count(sample11,minOverlap=3,summits=200)

Re-centering peaks...
Error in heights * sapply(called, function(x) x)

-In the second option of analysis:

 sample1<-dba.count(sample1,peaks=sample1$masks$Consensus,summits=200)

Re-centering peaks...
Error in 1:ncol(called) : argument of length 0

-in the third option:

sample1<-dba.count(sample1,peaks=sample1$masks$Consensus)

Error in pv$binding[, colnum]

Thank you in advance for your help

 

Paola

 

ADD COMMENT
0
Entering edit mode
Rory Stark ★ 5.1k
@rory-stark-5741
Last seen 8 days ago
Cambridge, UK

Hmmn. Did you include your original line of code adding a new sample? It won't work with that line in there!

Specifically, I expect the following to work:

> sample1 <- dba(sampleSheet="sample1.csv") #or whatever the samplesheet is called
> sample1 <- dba.count(sample1, minOverlap=3, summits=200)

These two lines work on my test data.

If it still doesn't work, try the following:

> sample1 <- dba.count(sample1, minOverlap=3, summits=TRUE)
> sample1 <- dba.count(sample1, peaks=NULL, summits=200)

If this still doesn't work, email me a link to the sample1 object after the first call to dba.count() with summits=TRUE, and I can track it down further. Please also email me the output from sessionInfo() so I can check versions of other packages.

Cheers-

Rory

ADD COMMENT
0
Entering edit mode
pabla86 • 0
@pabla86-12186
Last seen 7.2 years ago

Hi,

I've run the R code as you suggested, obtaining the following error message:

>sample1 <- dba.count(sample1, minOverlap=3, summits=TRUE)

Sample: IP_mapped_remove_dup_bam.bam125
Sample: input_mapped_remove_dup_bam.bam125
Error in pv$binding[, colnum] : incorrect number of dimensions
 Warning messages:
1: In dba.multicore.init(DBA$config) :
  Parallel execution unavailable: executing serially.
2: Unable to TMM normalize -- not enough peaksets

 

I'm sending you an email with the csv file.

Thank you in advance for your help.

Paola

 

ADD COMMENT

Login before adding your answer.

Traffic: 597 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6