Search
Question: Gviz Alignment Tracks with non-UCSC reference chromosome names
0
gravatar for laura.k.hilton
12 weeks ago by
laura.k.hilton0 wrote:

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
ADD COMMENTlink modified 12 weeks ago • written 12 weeks ago by laura.k.hilton0

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 REPLYlink written 12 weeks ago by florian.hahne@novartis.com1.6k

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 REPLYlink modified 12 weeks ago • written 12 weeks ago by laura.k.hilton0

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 REPLYlink written 12 weeks ago by laura.k.hilton0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 272 users visited in the last hour