Dear all,

I have an issue generating consensus peaks using the DiffBind package. I generate consensus peaks using the dba.peakset() function

dbaObj <-  dba(sampleSheet = "samples.csv")  

data.peakset <- dba.peakset(dbaObj, consensus = DBA_CONDITION, minOverlap = 2)

This generates columns for the consensus peaks in

data.peakset$allvectors

However, when I try to retrieve the peaks using dba.peakset with bRetrieve =TRUE.

data.peakset <- dba.peakset(dbaObj, consensus = DBA_CONDITION,
                            minOverlap = 2, bRetrieve = TRUE)

there are no consensus columns in the output. What is the supported way to get them?

Best wishes,

Georg

> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] grid      stats4    parallel  stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] DiffBind_1.12.1         limma_3.22.1            ChIPpeakAnno_2.16.2    
 [4] AnnotationDbi_1.28.1    Biobase_2.26.0          RSQLite_1.0.0          
 [7] DBI_0.3.1               biomaRt_2.22.0          VennDiagram_1.6.9      
[10] chipseq_1.16.0          ShortRead_1.24.0        GenomicAlignments_1.2.1
[13] Rsamtools_1.18.2        BiocParallel_1.0.0      BSgenome_1.34.0        
[16] rtracklayer_1.26.2      Biostrings_2.34.0       XVector_0.6.0          
[19] GenomicRanges_1.18.3    GenomeInfoDb_1.2.3      IRanges_2.0.0          
[22] S4Vectors_0.4.0         BiocGenerics_0.12.1    

loaded via a namespace (and not attached):
 [1] amap_0.8-12            base64enc_0.1-2        BatchJobs_1.5         
 [4] BBmisc_1.8             bitops_1.0-6           brew_1.0-6            
 [7] caTools_1.17.1         checkmate_1.5.0        codetools_0.2-9       
[10] compiler_3.1.1         digest_0.6.4           edgeR_3.8.5           
[13] fail_1.2               foreach_1.4.2          gdata_2.13.3          
[16] GenomicFeatures_1.18.2 GO.db_3.0.0            gplots_2.14.2         
[19] gtools_3.4.1           hwriter_1.3.2          iterators_1.0.7       
[22] KernSmooth_2.23-13     lattice_0.20-29        latticeExtra_0.6-26   
[25] MASS_7.3-35            multtest_2.22.0        RColorBrewer_1.0-5    
[28] RCurl_1.95-4.4         sendmailR_1.2-1        splines_3.1.1         
[31] stringr_0.6.2          survival_2.37-7        tools_3.1.1           
[34] XML_3.98-1.1           zlibbioc_1.12.0       

Hi Georg-

I can see how the interface to dba.peakset() is confusing in this case. The first time you call it, you are using it to add consensus peaksets to an existing DBA object, and it returns a new DBA object with the consensus peaksets added. The second time you are calling it, with bRetrieve=TRUE, you are retrieving a peakset, not a DBA object. Nothing is added, and the consensus and minOverlap parameters are ignored. It is not entirely clear from the documentation that the other parameters are ignored whenever bRetrieve=TRUE..

If you want to retrieve a peakset that includes metadata columns for the consensus peaksets, all you need to do is call dba.peakset twice in a row:

> consensusObj <- dba.peakset(dbaObj, consensus=DBA_CONDITION, minOverlap=2)
> data.peakset <- dba.peakset(consensusObj, bRetrieve=TRUE)

This will return a GRanges object with all the peaks, and metadata columns for all the original peaksets plus any consensus peaksets added by the first call to dba.peakset().

If you only want metadata columns for the consensus peaksets without the original peaksets as well, the sequence is:

> consensusObj <- dba.peakset(dbaObj, consensus=DBA_CONDITION, minOverlap=2)
> consensusObj <- dba(consensusObj, mask=consensusObj$masks$Consensus)
> data.peakset <- dba.peakset(consensusObj, bRetrieve=TRUE)

I'll have a look at changing the behavior of dba.peakset() so it honors the other parameters even when bRetrieve=TRUE, returning a peakset based on an internal, "hidden" DBA object.

Cheers-

Rory