Error in csaw::windowCounts when using readParam
Entering edit mode
Last seen 12 months ago

I'm getting an error from csaw::windowCounts when I use a readParam object in the params argument This does not change when I explicitly load the library(GenomeInfoDb) which exports seqinfo. Is this a bug or am I missing something?

This is the code producing the error:

blacklist <- rtracklayer::import.bed(
rpar <- csaw::readParam(
    minq = 20, discard = blacklist, restrict = c(as.character(as.roman(1:6)), "X")

# reproducible example Bam file
h3k27me3data <- H3K27me3Data()
bam <- h3k27me3data$Path[1]
# bam <- "out/alignments/ABC_123_sorted.bam" # my local bam

csaw::windowCounts(bam, param = rpar)
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘seqinfo’ for signature ‘"NULL"’

just calling seqinfo on blacklist works fine:

Seqinfo object with 6 sequences from an unspecified genome; no seqlengths:
  seqnames seqlengths isCircular genome
  chrIII           NA         NA   <NA>
  chrII            NA         NA   <NA>
  chrI             NA         NA   <NA>
  chrIV            NA         NA   <NA>
  chrV             NA         NA   <NA>
  chrX             NA         NA   <NA>


R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/
LAPACK: /usr/lib/x86_64-linux-gnu/

 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               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    LC_PAPER=en_GB.UTF-8       LC_NAME=C                 

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

other attached packages:
[1] testthat_2.3.2 devtools_2.3.1 usethis_1.6.1  reprex_0.3.0  

loaded via a namespace (and not attached):
 [1] Biobase_2.46.0              httr_1.4.2                  edgeR_3.28.1                pkgload_1.1.0              
 [5] bit64_4.0.2                 assertthat_0.2.1            askpass_1.1                 stats4_3.6.3               
 [9] BiocFileCache_1.10.2        blob_1.2.1                  GenomeInfoDbData_1.2.2      Rsamtools_2.2.3            
[13] remotes_2.2.0               progress_1.2.2              sessioninfo_1.1.1           pillar_1.4.6               
[17] RSQLite_2.2.0               backports_1.1.9             lattice_0.20-38             limma_3.42.2               
[21] glue_1.4.2                  digest_0.6.25               GenomicRanges_1.38.0        XVector_0.26.0             
[25] Matrix_1.2-18               XML_3.99-0                  pkgconfig_2.0.3             csaw_1.20.0                
[29] biomaRt_2.42.1              zlibbioc_1.32.0             purrr_0.3.4                 processx_3.4.3             
[33] BiocParallel_1.20.1         tibble_3.0.3                openssl_1.4.2               generics_0.0.2             
[37] IRanges_2.20.2              ellipsis_0.3.1              withr_2.2.0                 SummarizedExperiment_1.16.1
[41] GenomicFeatures_1.38.2      BiocGenerics_0.32.0         cli_2.0.2                   magrittr_1.5               
[45] crayon_1.3.4                memoise_1.1.0               ps_1.3.4                    fs_1.5.0                   
[49] fansi_0.4.1                 pkgbuild_1.1.0              tools_3.6.3                 prettyunits_1.1.1          
[53] hms_0.5.3                   lifecycle_0.2.0             matrixStats_0.56.0          stringr_1.4.0              
[57] S4Vectors_0.24.4            locfit_1.5-9.4              DelayedArray_0.12.3         AnnotationDbi_1.48.0       
[61] callr_3.4.3                 Biostrings_2.54.0           compiler_3.6.3              GenomeInfoDb_1.22.1        
[65] rlang_0.4.7                 grid_3.6.3                  RCurl_1.98-1.2              rstudioapi_0.11            
[69] rappdirs_0.3.1              bitops_1.0-6                curl_4.3                    DBI_1.1.0                  
[73] R6_2.4.1                    GenomicAlignments_1.22.1    dplyr_1.0.2                 rtracklayer_1.46.0         
[77] bit_4.0.4                   rprojroot_1.3-2             desc_1.2.0                  stringi_1.4.6              
[81] parallel_3.6.3              Rcpp_1.0.5                  vctrs_0.3.3                 tidyselect_1.1.0           
[85] dbplyr_1.4.4
csaw GenomeInfoDb • 298 views
Entering edit mode
Aaron Lun ★ 27k
Last seen 14 hours ago
The city by the bay

The problem is that restrict= is set to roman numbers (presumably for Drosophila?) but obviously the data itself is for mouse, which makes the function very confused. Arguably the failure could be more graceful, but regardless, the function won't be doing anything sensible with these settings.

Entering edit mode

Ah, of course, the chromosome names don't match (it's C. elegans data) I forgot to prefix with 'chr'. paste0("chr", c(as.character(as.roman(1:6)), "X")) works fine. The error message had me looking for namespace issues with the seqinfo function not being exported the child process this runs on in my actual code - but of course it's something simple! I imagine that mismatching chromosome names is a common enough issue that it would be nice to have an error message that checks for this though. Thanks


Login before adding your answer.

Traffic: 267 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6