Error in csaw::windowCounts when using readParam
1
1
Entering edit mode
@richardjacton-12268
Last seen 4.2 years 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(
    "https://github.com/Boyle-Lab/Blacklist/raw/master/lists/ce11-blacklist.v2.bed.gz"
)
rpar <- csaw::readParam(
    minq = 20, discard = blacklist, restrict = c(as.character(as.roman(1:6)), "X")
)

# reproducible example Bam file
library(chipseqDBData)
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(blacklist)
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>

sessionInfo:

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/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [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                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=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 • 1.4k views
ADD COMMENT
0
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 4 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.

ADD COMMENT
0
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

ADD REPLY

Login before adding your answer.

Traffic: 806 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