Question

DiffBind dba.count error pv$merge(bed): attempt to apply non-function

0

Entering edit mode

chadnewton • 0

@chadnewton-20911

Last seen 6.0 years ago

I am attempting to perform differential peak analysis from ATAC-seq data using DiffBind package and am running into problems with the dba.count function.

I have 2 conditions (rapid vs slow) with 3 samples each for a total of 6 samples. I want to create a consensus peakset for both rapid and slow condition. I would like to identify punctate peaks, so I set the summits = 200.

cn <- dba(sampleSheet = "cn_atac.csv", peakFormat = "bed", minOverlap = 1, 
      config = data.frame(RunParallel=TRUE, reportInit='DBA', DataType=DBA_DATA_GRANGES,
                          minQCth=15, fragmentSize=125,
                          bCorPlot=FALSE, th=0.05, bUsePval=FALSE), 
      bRemoveM = TRUE, bRemoveRandom = TRUE, filter = 20)
cn.consensus = dba.peakset(cn, consensus = DBA_CONDITION, minOverlap = 1)
peaks = dba.peakset(cn.consensus, cn.consensus$masks$Consensus, bRetrieve = TRUE)
cn.cons.count = dba.count(cn, peaks = peaks, 
                      summits = 200, bRemoveDuplicates = FALSE, bScaleControl = TRUE,
                      bParallel = FALSE)

generates error:

 Sample: atac_files/c1.bam125 
    [W::bam_hdr_read] bgzf_check_EOF: No error
    [W::bam_hdr_read] bgzf_check_EOF: No error
    Sample: atac_files/c2.bam125 
    Sample: atac_files/c3.bam125 
    [W::bam_hdr_read] bgzf_check_EOF: No error
    [W::bam_hdr_read] bgzf_check_EOF: No error
    Sample: atac_files/c4.bam125 
    Sample: atac_files/c5.bam125 
    [W::bam_hdr_read] bgzf_check_EOF: No error
    [W::bam_hdr_read] bgzf_check_EOF: No error
    Sample: atac_files/c6.bam125 
    Re-centering peaks...
    Error in pv$merge(bed) : attempt to apply non-function
    In addition: Warning message:
      In dba.multicore.init(DBA$config) :
      Parallel execution unavailable: executing serially.

> traceback()
3: pv.Recenter(pv, summits, (numpeaks - numAdded + 1):numpeaks, 
       called)
2: pv.counts(DBA, peaks = peaks, minOverlap = minOverlap, defaultScore = score, 
       bLog = bLog, insertLength = fragmentSize, bOnlyCounts = T, 
       bCalledMasks = TRUE, minMaxval = filter, bParallel = bParallel, 
       bUseLast = bUseLast, bWithoutDupes = bRemoveDuplicates, bScaleControl = bScaleControl, 
       filterFun = filterFun, bLowMem = bUseSummarizeOverlaps, readFormat = readFormat, 
       summits = summits, minMappingQuality = mapQCth)
1: dba.count(cn, peaks = peaks, summits = 200, bRemoveDuplicates = FALSE, 
       bScaleControl = TRUE, bParallel = FALSE)

.

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

other attached packages:
 [1] DiffBind_2.12.0             SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
 [4] BiocParallel_1.17.18        matrixStats_0.54.0          Biobase_2.44.0             
 [7] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0         IRanges_2.18.0             
[10] S4Vectors_0.22.0            BiocGenerics_0.30.0        

loaded via a namespace (and not attached):
  [1] amap_0.8-17              colorspace_1.4-1         rjson_0.2.20             hwriter_1.3.2           
  [5] htmlTable_1.13.1         XVector_0.24.0           base64enc_0.1-3          rstudioapi_0.10         
  [9] ggrepel_0.8.1            bit64_0.9-7              AnnotationDbi_1.46.0     splines_3.6.0           
 [13] geneplotter_1.62.0       knitr_1.23               Formula_1.2-3            Rsamtools_2.0.0         
 [17] annotate_1.62.0          cluster_2.0.8            GO.db_3.8.2              pheatmap_1.0.12         
 [21] graph_1.62.0             BiocManager_1.30.4       compiler_3.6.0           httr_1.4.0              
 [25] GOstats_2.50.0           backports_1.1.4          assertthat_0.2.1         Matrix_1.2-17           
 [29] lazyeval_0.2.2           limma_3.40.2             htmltools_0.3.6          acepack_1.4.1           
 [33] prettyunits_1.0.2        tools_3.6.0              gtable_0.3.0             glue_1.3.1              
 [37] GenomeInfoDbData_1.2.1   Category_2.50.0          systemPipeR_1.18.0       dplyr_0.8.1             
 [41] batchtools_0.9.11        rappdirs_0.3.1           ShortRead_1.42.0         Rcpp_1.0.1              
 [45] Biostrings_2.52.0        gdata_2.18.0             rtracklayer_1.44.0       xfun_0.7                
 [49] stringr_1.4.0            gtools_3.8.1             XML_3.98-1.19            edgeR_3.26.3            
 [53] zlibbioc_1.30.0          scales_1.0.0             BSgenome_1.52.0          VariantAnnotation_1.30.1
 [57] hms_0.4.2                RBGL_1.60.0              RColorBrewer_1.1-2       yaml_2.2.0              
 [61] gridExtra_2.3            memoise_1.1.0            ggplot2_3.1.1            biomaRt_2.40.0          
 [65] rpart_4.1-15             latticeExtra_0.6-28      stringi_1.4.3            RSQLite_2.1.1           
 [69] genefilter_1.66.0        checkmate_1.9.3          GenomicFeatures_1.36.0   caTools_1.17.1.2        
 [73] rlang_0.3.4              pkgconfig_2.0.2          bitops_1.0-6             lattice_0.20-38         
 [77] purrr_0.3.2              labeling_0.3             htmlwidgets_1.3          GenomicAlignments_1.20.0
 [81] bit_1.1-14               tidyselect_0.2.5         GSEABase_1.46.0          AnnotationForge_1.26.0  
 [85] plyr_1.8.4               magrittr_1.5             DESeq2_1.24.0            R6_2.4.0                
 [89] gplots_3.0.1.1           Hmisc_4.2-0              base64url_1.4            DBI_1.0.0               
 [93] pillar_1.4.0             foreign_0.8-71           withr_2.1.2              survival_2.44-1.1       
 [97] RCurl_1.95-4.12          nnet_7.3-12              tibble_2.1.1             crayon_1.3.4            
[101] KernSmooth_2.23-15       progress_1.2.2           locfit_1.5-9.1           grid_3.6.0              
[105] data.table_1.12.2        blob_1.1.1               Rgraphviz_2.28.0         digest_0.6.19           
[109] xtable_1.8-4             brew_1.0-6               munsell_0.5.0

Much appreciated! Chad

DiffBind • 1.6k views

ADD COMMENT • link updated 6.2 years ago by Rory Stark ★ 5.2k • written 6.4 years ago by chadnewton • 0

0

Entering edit mode

Did you figure this out? I'm having much the same issue. There's no discernible difference in our code so I wont dump that on here, but I have included my session info below if anyone else needs it.

Another similarity is that I am trying to use this for ATACseq. With that I do not having any inputs. Could this contribute?

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

Random number generation:
 RNG:     Mersenne-Twister 
 Normal:  Inversion 
 Sample:  Rounding 

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] DiffBind_2.12.0             SummarizedExperiment_1.14.0 DelayedArray_0.10.0         BiocParallel_1.18.0        
 [5] matrixStats_0.54.0          Biobase_2.44.0              GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
 [9] IRanges_2.18.1              S4Vectors_0.22.0            BiocGenerics_0.30.0         BiFET_1.4.0                

loaded via a namespace (and not attached):
 [1] Category_2.50.0          bitops_1.0-6             bit64_0.9-7              RColorBrewer_1.1-2       progress_1.2.2          
 [6] httr_1.4.0               Rgraphviz_2.28.0         backports_1.1.4          tools_3.6.0              R6_2.4.0                
[11] KernSmooth_2.23-15       DBI_1.0.0                lazyeval_0.2.2           colorspace_1.4-1         withr_2.1.2             
[16] tidyselect_0.2.5         prettyunits_1.0.2        bit_1.1-14               compiler_3.6.0           graph_1.62.0            
[21] rtracklayer_1.44.0       checkmate_1.9.3          caTools_1.17.1.2         scales_1.0.0             genefilter_1.66.0       
[26] RBGL_1.60.0              rappdirs_0.3.1           stringr_1.4.0            digest_0.6.19            Rsamtools_2.0.0         
[31] AnnotationForge_1.26.0   XVector_0.24.0           pkgconfig_2.0.2          BSgenome_1.52.0          limma_3.40.2            
[36] rlang_0.3.4              rstudioapi_0.10          RSQLite_2.1.1            GOstats_2.50.0           hwriter_1.3.2           
[41] gtools_3.8.1             dplyr_0.8.1              VariantAnnotation_1.30.1 RCurl_1.95-4.12          magrittr_1.5            
[46] GO.db_3.8.2              GenomeInfoDbData_1.2.1   Matrix_1.2-17            Rcpp_1.0.1               munsell_0.5.0           
[51] stringi_1.4.3            yaml_2.2.0               edgeR_3.26.4             zlibbioc_1.30.0          gplots_3.0.1.1          
[56] plyr_1.8.4               grid_3.6.0               blob_1.1.1               ggrepel_0.8.1            gdata_2.18.0            
[61] crayon_1.3.4             lattice_0.20-38          splines_3.6.0            Biostrings_2.52.0        GenomicFeatures_1.36.1  
[66] annotate_1.62.0          poibin_1.3               hms_0.4.2                batchtools_0.9.11        locfit_1.5-9.1          
[71] knitr_1.23               pillar_1.4.1             rjson_0.2.20             systemPipeR_1.18.1       base64url_1.4           
[76] biomaRt_2.40.0           XML_3.98-1.20            glue_1.3.1               ShortRead_1.42.0         latticeExtra_0.6-28     
[81] data.table_1.12.2        gtable_0.3.0             purrr_0.3.2              amap_0.8-17              assertthat_0.2.1        
[86] ggplot2_3.1.1            xfun_0.7                 xtable_1.8-4             survival_2.44-1.1        pheatmap_1.0.12         
[91] tibble_2.1.3             GenomicAlignments_1.20.0 AnnotationDbi_1.46.0     memoise_1.1.0            brew_1.0-6              
[96] GSEABase_1.46.0

ADD REPLY • link 6.4 years ago matthew.paul.2006 ▴ 10

score 0 · Answer 1 · 2019-08-22

I'm suspicous that some of the files generate a [W::bam_hdr_read] bgzf_check_EOF: No error warning and some do not.

Can you try calling dba.count() a) without setting summits or b) setting summits=TRUE? Does it work in either of these cases? If it works with summits=TRUE, you can run dba.count() again with peaks=NULL and summits=200, and see if that works.

If it works with summits=TRUE, but fails when you call it again with peaks=NULL and summits=200, please send me a link to the DBA object after calling with with summits=TRUE so I can have a look.

-Rory