Question: DiffBind dba.peakset() and consensus peaks
0
gravatar for Georg Otto
5.0 years ago by
Georg Otto120
United Kingdom
Georg Otto120 wrote:

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       

 

chipseq diffbind • 2.7k views
ADD COMMENTlink modified 5.0 years ago by Rory Stark3.0k • written 5.0 years ago by Georg Otto120
Answer: DiffBind dba.peakset() and consensus peaks
3
gravatar for Rory Stark
5.0 years ago by
Rory Stark3.0k
CRUK, Cambridge, UK
Rory Stark3.0k wrote:

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

 

ADD COMMENTlink written 5.0 years ago by Rory Stark3.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 220 users visited in the last hour