Gviz Alignment Tracks with non-UCSC reference chromosome names
0
0
Entering edit mode
@laurakhilton-17332
Last seen 5.6 years ago

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
gviz software error • 1.4k views
ADD COMMENT
0
Entering edit mode

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.

 

 

ADD REPLY
0
Entering edit mode

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: 

library(Gviz)

afrom <- 161645010
ato <- 161645990

# For bam file with chromosome="1" notation

bamFile <- "02-13135.normal.FCGR.GRCh37.bam"

options(ucscChromosomeNames=FALSE)
atrack <- AlignmentsTrack(bamFile, genome="hg19", chromosome="1", start=afrom, end=ato, type="coverage")

plotTracks(atrack, from=afrom, to=ato, chromosome="1")

# Produces error:
# Error in .normarg_at2(at, x) : 
#   some ranges in 'at' are off-limits with respect to their corresponding sequence in 'x'

# For bam file with chromosome="chr1" notation

bamFile <- "02-13135.normal.FCGR.GRCh37.chr.bam"

atrack <- AlignmentsTrack(bamFile, genome="hg19", chromosome="chr1", start=afrom, end=ato, type="coverage")

plotTracks(atrack, from=afrom, to=ato, chromosome="chr1")

#Error: Crashes R

And here is the session info from my most recent test of this script: 

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.7 (Final)

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

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

other attached packages:
[1] Gviz_1.18.2          GenomicRanges_1.26.4 GenomeInfoDb_1.10.3  IRanges_2.8.2        S4Vectors_0.12.2     BiocGenerics_0.20.0 

loaded via a namespace (and not attached):
 [1] base64_2.0                    Rcpp_0.12.11                  biovizBase_1.22.0             lattice_0.20-34              
 [5] Rsamtools_1.26.1              Biostrings_2.42.1             digest_0.6.12                 mime_0.5                     
 [9] R6_2.2.0                      plyr_1.8.4                    acepack_1.4.1                 RSQLite_1.1-1                
[13] BiocInstaller_1.24.0          httr_1.2.1                    ggplot2_2.2.1                 zlibbioc_1.20.0              
[17] rlang_0.1.1                   GenomicFeatures_1.26.0        lazyeval_0.2.0                data.table_1.10.0            
[21] rpart_4.1-10                  Matrix_1.2-7.1                splines_3.3.2                 BiocParallel_1.8.1           
[25] AnnotationHub_2.6.4           stringr_1.2.0                 foreign_0.8-67                RCurl_1.95-4.8               
[29] biomaRt_2.30.0                munsell_0.5.0                 shiny_0.14.2                  httpuv_1.3.3                 
[33] rtracklayer_1.34.1            htmltools_0.3.5               nnet_7.3-12                   openssl_0.9.5                
[37] SummarizedExperiment_1.4.0    tibble_1.3.3                  gridExtra_2.2.1               htmlTable_1.7                
[41] interactiveDisplayBase_1.12.0 Hmisc_4.0-1                   matrixStats_0.52.2            XML_3.98-1.5                 
[45] GenomicAlignments_1.10.0      bitops_1.0-6                  xtable_1.8-2                  gtable_0.2.0                 
[49] DBI_0.5-1                     magrittr_1.5                  scales_1.0.0                  stringi_1.1.5                
[53] XVector_0.14.0                latticeExtra_0.6-28           Formula_1.2-1                 RColorBrewer_1.1-2           
[57] ensembldb_1.6.2               tools_3.3.2                   dichromat_2.0-0               BSgenome_1.42.0              
[61] Biobase_2.34.0                yaml_2.1.14                   survival_2.39-5               AnnotationDbi_1.36.0         
[65] colorspace_1.3-2              cluster_2.0.5                 memoise_1.0.0                 VariantAnnotation_1.20.2     
[69] knitr_1.15.1         

 

ADD REPLY
0
Entering edit mode

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...

ADD REPLY

Login before adding your answer.

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