Hello!
I am trying to plotTracks
reads mapped to a ChIP-seq peak region for a set of 6 BAM files (3 ChIP, 3 inputs) that contain paired-end reads. If I subset to certain samples, it works as per this screenshot: https://www.dropbox.com/s/s71s9wu8v8vpshb/Screenshot%202018-06-19%2014.09.19.png?dl=1
However, certain samples (1 ChIP, 2 inputs) cause the error:
Error in .normarg_at2(at, x) : some ranges in 'at' are off-limits with respect to their corresponding sequence in 'x'
It's a bit challenging to share a minimal reproducible example in the post, so I have linked to my public Dropbox for now.
https://www.dropbox.com/s/htpmsqfadut3iyr/Gviz.zip?dl=1
This link points to a ZIP archive (13 KB) that contains small files to reproduce my issue:
- region.bed: the peak region
- subset.bam and the accompanying *bai index: reads that overlap the region defined above. I generated those files using
bedtools intersect
to produce a small file that I can share to reproduce the issue.- Note that this file produces the exact same error as the full-size file, even though it probable misses some reads outside the region that are paired with a mate inside the region.
- issue.R: an R script that uses the files above and produces the error.
For the record, issue.R looks like this:
require(rtracklayer) topRegion <- import.bed("region.bed") require(Rsamtools) sbp <- ScanBamParam( which=GRanges( seqnames = seqnames(topRegion), ranges=IRanges( start=start(topRegion)-width(topRegion), end = end(topRegion)+width(topRegion)) ) ) require(Gviz) testTrack <- AlignmentsTrack("subset.bam", isPaired=TRUE, name = "Test") plotTracks( list(testTrack), from=start(topRegion), to=end(topRegion), chromosome=seqnames(topRegion), extend.left = width(topRegion), extend.right = width(topRegion))
Here is the stack trace that I can see in RStudio:
Error in .normarg_at2(at, x) : some ranges in 'at' are off-limits with respect to their corresponding sequence in 'x' 21. stop(.wrap_msg("some ranges in 'at' are off-limits with respect to ", "their corresponding sequence in 'x'")) 20. .normarg_at2(at, x) 19. replaceAt(x, at, value = value) 18. replaceAt(x, at, value = value) 17. sequenceLayer(reads$seq, reads$cigar) 16. x@stream(x@reference, subRegion) 15. mcols(data) 14. colnames(mcols(data)) 13. eval(quote(list(...)), env) 12. eval(quote(list(...)), env) 11. eval(quote(list(...)), env) 10. standardGeneric("paste") 9. paste(colnames(mcols(data)), "orig", sep = "__.__") 8. .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping) 7. .local(x, ...) 6. subset(x, from = from, to = to, chromosome = chromosome, sort = FALSE, stacks = FALSE, use.defaults = FALSE) 5. subset(x, from = from, to = to, chromosome = chromosome, sort = FALSE, stacks = FALSE, use.defaults = FALSE) 4. FUN(X[[i]], ...) 3. lapply(trackList, function(x) { chromosome(x) <- chromosome subset(x, from = from, to = to, chromosome = chromosome, sort = FALSE, stacks = FALSE, use.defaults = FALSE) ... 2. lapply(trackList, function(x) { chromosome(x) <- chromosome subset(x, from = from, to = to, chromosome = chromosome, sort = FALSE, stacks = FALSE, use.defaults = FALSE) ... 1. plotTracks(list(testTrack), from = start(topRegion), to = end(topRegion), chromosome = seqnames(topRegion), extend.left = width(topRegion), extend.right = width(topRegion))
Here is my session information:
> sessionInfo() R version 3.5.0 (2018-04-23) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS High Sierra 10.13.5 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib locale: [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 attached base packages: [1] grid parallel stats4 stats graphics grDevices utils datasets methods base other attached packages: [1] Gviz_1.24.0 Rsamtools_1.32.0 Biostrings_2.48.0 XVector_0.20.0 rtracklayer_1.40.3 [6] GenomicRanges_1.32.3 GenomeInfoDb_1.16.0 IRanges_2.14.10 S4Vectors_0.18.3 BiocGenerics_0.26.0 loaded via a namespace (and not attached): [1] Biobase_2.40.0 httr_1.3.1 bit64_0.9-7 splines_3.5.0 [5] Formula_1.2-3 assertthat_0.2.0 latticeExtra_0.6-28 blob_1.1.1 [9] BSgenome_1.48.0 GenomeInfoDbData_1.1.0 yaml_2.1.19 progress_1.1.2 [13] pillar_1.2.3 RSQLite_2.1.1 backports_1.1.2 lattice_0.20-35 [17] biovizBase_1.28.0 digest_0.6.15 RColorBrewer_1.1-2 checkmate_1.8.5 [21] colorspace_1.3-2 htmltools_0.3.6 Matrix_1.2-14 plyr_1.8.4 [25] XML_3.98-1.11 biomaRt_2.36.1 zlibbioc_1.26.0 scales_0.5.0 [29] BiocParallel_1.14.1 htmlTable_1.12 tibble_1.4.2 AnnotationFilter_1.4.0 [33] ggplot2_2.2.1 SummarizedExperiment_1.10.1 GenomicFeatures_1.32.0 nnet_7.3-12 [37] lazyeval_0.2.1 survival_2.42-3 magrittr_1.5 memoise_1.1.0 [41] foreign_0.8-70 tools_3.5.0 data.table_1.11.4 prettyunits_1.0.2 [45] matrixStats_0.53.1 stringr_1.3.1 munsell_0.4.3 cluster_2.0.7-1 [49] DelayedArray_0.6.0 AnnotationDbi_1.42.1 ensembldb_2.4.1 compiler_3.5.0 [53] rlang_0.2.1 RCurl_1.95-4.10 dichromat_2.0-0 rstudioapi_0.7 [57] VariantAnnotation_1.26.0 htmlwidgets_1.2 bitops_1.0-6 base64enc_0.1-3 [61] gtable_0.2.0 curl_3.2 DBI_1.0.0 R6_2.2.2 [65] GenomicAlignments_1.16.0 gridExtra_2.3 knitr_1.20 bit_1.1-14 [69] Hmisc_4.1-1 ProtGenerics_1.12.0 stringi_1.2.2 Rcpp_0.12.17 [73] rpart_4.1-13 acepack_1.4.1
I just realized I just posted almost exactly the same question here: Gviz Alignment Tracks with non-UCSC reference chromosome names