I'm having trouble plotting an alignment track from sequence data that was aligned to hg19 with non-UCSC chromosome names (i.e. "1", not "chr1"). I suspect the problem is to do with the chromosome names, but I am not sure. My code is as follows:
afrom <- 161564010 ato <- 161645990
bamFile <- "/path/to/bamfile"
options(ucscChromosomeNames = FALSE) atrack <- AlignmentsTrack(bamFile, isPaired=TRUE, genome = "hg19")
plotTracks(atrack, from=afrom, to=ato, chromosome="1")
I get the following error after executing the plotTracks function:
Error in .normarg_at2(at, x) : some ranges in 'at' are off-limits with respect to their corresponding sequence in 'x'
> traceback() 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(atrack, from = afrom, to = ato, chromosome = "1")
If instead I run it using plotTracks(atrack, from=afrom, to=ato, chromosome="chr1")
, the error vanishes and it outputs a plot that is blank except for the plot name on the Y-axis, obviously because my bam file doesn't contain any sequences on "chr1."
Any suggestions about how I can resolve this? I can send my bam file for debugging if needed. Thanks.
sessionInfo() R version 3.4.1 (2017-06-30) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.2 LTS Matrix products: default BLAS: /usr/lib/libblas/libblas.so.3.6.0 LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 locale: [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8 LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8 [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid parallel stats4 stats graphics grDevices utils datasets methods base other attached packages: [1] readr_1.1.1 biomaRt_2.34.2 Homo.sapiens_1.3.1 org.Hs.eg.db_3.5.0 [5] GO.db_3.5.0 OrganismDbi_1.20.0 TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.30.3 [9] AnnotationDbi_1.40.0 VariantAnnotation_1.24.5 Rsamtools_1.30.0 Biostrings_2.46.0 [13] XVector_0.18.0 SummarizedExperiment_1.8.1 DelayedArray_0.4.1 matrixStats_0.54.0 [17] Biobase_2.38.0 Gviz_1.22.3 GenomicRanges_1.30.3 GenomeInfoDb_1.14.0 [21] IRanges_2.12.0 S4Vectors_0.16.0 BiocGenerics_0.24.0 loaded via a namespace (and not attached): [1] ProtGenerics_1.10.0 bitops_1.0-6 bit64_0.9-7 RColorBrewer_1.1-2 progress_1.2.0 [6] httr_1.3.1 tools_3.4.1 backports_1.1.2 R6_2.2.2 rpart_4.1-11 [11] Hmisc_4.1-1 DBI_1.0.0 lazyeval_0.2.1 colorspace_1.3-2 nnet_7.3-12 [16] tidyselect_0.2.4 gridExtra_2.3 prettyunits_1.0.2 RMySQL_0.10.15 curl_3.2 [21] bit_1.1-14 compiler_3.4.1 graph_1.56.0 htmlTable_1.12 rtracklayer_1.38.3 [26] scales_1.0.0 checkmate_1.8.5 RBGL_1.54.0 stringr_1.3.1 digest_0.6.17 [31] foreign_0.8-69 base64enc_0.1-3 dichromat_2.0-0 pkgconfig_2.0.2 htmltools_0.3.6 [36] ensembldb_2.2.2 BSgenome_1.46.0 htmlwidgets_1.2 rlang_0.2.2 rstudioapi_0.7 [41] RSQLite_2.1.1 BiocInstaller_1.28.0 shiny_1.1.0 bindr_0.1.1 BiocParallel_1.12.0 [46] acepack_1.4.1 dplyr_0.7.6 RCurl_1.95-4.11 magrittr_1.5 GenomeInfoDbData_1.0.0 [51] Formula_1.2-3 Matrix_1.2-10 Rcpp_0.12.18 munsell_0.5.0 stringi_1.2.4 [56] yaml_2.2.0 zlibbioc_1.24.0 plyr_1.8.4 AnnotationHub_2.10.1 blob_1.1.1 [61] promises_1.0.1 crayon_1.3.4 lattice_0.20-35 splines_3.4.1 hms_0.4.2 [66] knitr_1.20 pillar_1.3.0 XML_3.98-1.16 glue_1.3.0 biovizBase_1.26.0 [71] latticeExtra_0.6-28 data.table_1.11.4 httpuv_1.4.5 gtable_0.2.0 purrr_0.2.5 [76] assertthat_0.2.0 ggplot2_3.0.0 mime_0.5 xtable_1.8-3 AnnotationFilter_1.2.0 [81] later_0.7.4 survival_2.41-3 tibble_1.4.2 GenomicAlignments_1.14.2 memoise_1.1.0 [86] bindrcpp_0.2.2 cluster_2.0.6 interactiveDisplayBase_1.16.0
Could you create a small reproducible example, maybe with a stripped down version of the BAM file?
Regardless of the chromosome issue, this should fail more gracefully. Could you please try and set
options(ucscChromosomeNames=FALSE)
That might help with the non-standard chromosome identifiers.
Hi Florian,
Thank you for the prompt reply. I did indeed have
options(ucscChromosomeNames=FALSE)
set when I was running this.I also tried generating a bam file that used UCSC "chr1" notation instead of Entrez "1" notation. This indeed removed the error, and instead R crashes whenever I try to run my script. It crashes regardless of whether I use Rstudio or run it at the command line. So it appears I've solved the error I mentioned before, but Gviz crashes whenever I try to plot my alignment track.
Here is a link to a Dropbox folder with my BAM files: https://www.dropbox.com/sh/acyd1mhp81u7i95/AABRD11csl434OZq-ejtlcAja?dl=0
Note that the bam file with ".chr" in the filename has the UCSC chromosome notation. The other has Entrez notation.
And here is my minimal reproducible R script:
And here is the session info from my most recent test of this script:
It seems to be an issue with my bam file. I ran it again on a whole genome bam file (the files I sent you were pared down to a small region). Gviz worked just fine in this instance.
No idea what could possibly be wrong with my bam file though...