Can you compute summary stats of coverage directly when importing using rtracklayer?
0
0
Entering edit mode
David • 0
@david-24104
Last seen 2.7 years ago
United Kingdom

I am trying to obtain coverage summary stats (e.g. mean/sum) from a BigWig across a set of ranges. I would like to know if I'm currently doing this in the fastest possible way using rtracklayer::import.

If I set as = "GRanges", this gives one score per range in the BigWig (rather than the desired 1 score per input range).

If I set as = "NumericList", I can obtain the coverage across individual bases then compute the summary stats for each range using e.g. BiocGenerics::mean. However, I was wondering if there exists a (faster) method of calculating this directly using rtracklayer::import? I can't seem to find anything in the vignette or function documentation.

test_path <- system.file("tests", package = "rtracklayer")
test_bw <- file.path(test_path, "test.bw")

ranges <- GenomicRanges::GRanges(c("chr2:1-400", "chr2:500-1000"))

gr <- rtracklayer::import(test_bw, which = ranges, as = "GRanges")
gr 
#> GRanges object with 5 ranges and 1 metadata column:
#>       seqnames    ranges strand |     score
#>          <Rle> <IRanges>  <Rle> | <numeric>
#>   [1]     chr2     1-300      * |     -1.00
#>   [2]     chr2   301-400      * |     -0.75
#>   [3]     chr2   500-600      * |     -0.75
#>   [4]     chr2   601-900      * |     -0.50
#>   [5]     chr2  901-1000      * |     -0.25
#>   -------
#>   seqinfo: 2 sequences from an unspecified genome

nl <- rtracklayer::import(test_bw, which = ranges, as = "NumericList")
BiocGenerics::mean(nl)
#> [1] -0.937500 -0.500499

sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/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] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] pillar_1.4.6                compiler_4.0.3             
#>  [3] highr_0.8                   GenomeInfoDb_1.26.0        
#>  [5] XVector_0.30.0              MatrixGenerics_1.2.0       
#>  [7] bitops_1.0-6                tools_4.0.3                
#>  [9] zlibbioc_1.36.0             digest_0.6.27              
#> [11] lattice_0.20-41             evaluate_0.14              
#> [13] lifecycle_0.2.0             tibble_3.0.4               
#> [15] pkgconfig_2.0.3             rlang_0.4.8                
#> [17] Matrix_1.2-18               reprex_0.3.0.9001          
#> [19] DelayedArray_0.16.0         rstudioapi_0.11            
#> [21] yaml_2.2.1                  parallel_4.0.3             
#> [23] xfun_0.19                   GenomeInfoDbData_1.2.4     
#> [25] rtracklayer_1.50.0          styler_1.3.2               
#> [27] stringr_1.4.0               knitr_1.30                 
#> [29] Biostrings_2.58.0           fs_1.5.0                   
#> [31] vctrs_0.3.4                 S4Vectors_0.28.0           
#> [33] IRanges_2.24.0              grid_4.0.3                 
#> [35] stats4_4.0.3                Biobase_2.50.0             
#> [37] XML_3.99-0.5                BiocParallel_1.24.0        
#> [39] rmarkdown_2.5               purrr_0.3.4                
#> [41] magrittr_1.5                matrixStats_0.57.0         
#> [43] GenomicAlignments_1.26.0    backports_1.2.0            
#> [45] Rsamtools_2.6.0             ellipsis_0.3.1             
#> [47] htmltools_0.5.0             BiocGenerics_0.36.0        
#> [49] GenomicRanges_1.42.0        SummarizedExperiment_1.20.0
#> [51] stringi_1.5.3               RCurl_1.98-1.2             
#> [53] crayon_1.3.4

Created on 2020-11-09 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0.9001)

Coverage rtracklayer • 960 views
ADD COMMENT

Login before adding your answer.

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